#1. Data wrangling ##1.1 Loading the required packages
#install.packages("pacman")
pacman::p_load(tidyverse, ggplot2, brms, cmdstanr, broom, broom.mixed, tidybayes, ggthemes, emmeans, ggdist, gridExtra, bayesplot, see, patchwork, ggtext, extraDistr, gghalves, ggpp, rstan)
#remotes::install_github("stan-dev/loo")
# Setting colours for vidsalisations
colour1 <- "#97A275"
colour2 <- "#F2B23A"
colour3 <- "#F2711B"
colour4 <- "#8CBEB2"
colour5 <- "#E7C27D"
colour6 <- "#F2A775"
colour11 <- "#666E4F"
colour12 <- "#A38346"
colour13 <- "#9C4811"
##1.2 Unblinding the data
#assigning the directory so that the code can be universally read on any computer as long as the entire project file is shared
directory <- './2-data/'
#reading in the coded raw BORIS output
BORIS_output <- read.csv("./2-data/BORIS_output.csv")
#reading in the decoder
Video_codes <- read.csv("./2-data/Video_codes.csv")
#reading in shrimp species information
shrimp_species <- read.csv("./2-data/species.csv")
#binding the decoder to the BORIS output and revealing the treatment and stage
BORIS_output$code <- str_sub(BORIS_output$video_id, 2,)
BORIS_output$code <- as.factor(BORIS_output$code)
shrimp_species$code <- as.factor(shrimp_species$code)
Video_codes$code <- as.factor(Video_codes$code)
decoded_dataframe <- inner_join(inner_join(BORIS_output, Video_codes, by = c("code")), shrimp_species, by = c("code"))
# Find common codes across all data sets
common_codes <- Reduce(intersect, list(BORIS_output$code, shrimp_species$code, Video_codes$code))
# Combine all codes from the three data sets
all_codes <- unique(c(BORIS_output$code, shrimp_species$code, Video_codes$code))
# Find codes not common across all data sets
codes_not_common <- setdiff(all_codes, common_codes)
# Define a function to check which data sets contain each code
check_code_origin <- function(code) {
origins <- c()
if(code %in% BORIS_output$code) {
origins <- c(origins, "BORIS_output")
}
if(code %in% shrimp_species$code) {
origins <- c(origins, "shrimp_species")
}
if(code %in% Video_codes$code) {
origins <- c(origins, "Video_codes")
}
return(paste(origins, collapse=", "))
}
# Create a dataframe to store the results
results <- data.frame(code = character(0), origin = character(0), stringsAsFactors = FALSE)
# Loop through each code not common across all data sets
for(code in codes_not_common) {
origin <- check_code_origin(code)
results <- rbind(results, data.frame(code = as.character(code), origin = origin, stringsAsFactors = FALSE))
}
# Print the results
print(results)
## code origin
## 1 31 Video_codes
## 2 32 Video_codes
## 3 33 Video_codes
## 4 34 Video_codes
## 5 35 Video_codes
## 6 36 Video_codes
## 7 40 Video_codes
## 8 41 Video_codes
## 9 42 Video_codes
## 10 148 Video_codes
## 11 149 Video_codes
## 12 150 Video_codes
# Missing observations are identified exclusively within the 'Video_codes' dataframe. This situation suggests that all videos capturing both shrimp and gobies are consistently represented across the datasets. However, there were instances where burrows were recorded and assigned video codes but showed no activity from gobies or shrimp during these recordings. These instances led to the assumption that such burrows were inactive, resulting in their exclusion from the further analysis of the study.
##1.4 Cleaning the dataframe ###1.4.1 Column renaming
decoded_dataframe <- decoded_dataframe |>
rename(shrimp_species = species)
#renaming some of the variables and extracting individual variables from columns that have multiple variables contained in one string
decoded_dataframe <- decoded_dataframe |>
dplyr::rename(stage = video_section,
duration = total_duration_s,
behaviour = behavior) |>
dplyr::mutate(date = str_sub(date_file, , 6),
species = str_sub(subject, 5,),
duration = as.numeric(duration),
boat_id = str_sub(date_file, -1, -1),
video_num = as.numeric(str_sub(video_id, 2)),
burrow = ceiling(video_num/3),
individual = str_sub(subject, 1, 1),
stroke = factor(stroke),
boat_name = factor(boat_name)
)
#creating some unique identifiers to make grouping easier later on for reorganising of the dataframe
decoded_dataframe <- decoded_dataframe |>
unite("unique_id", 'burrow','individual', 'species', 'stage', remove = FALSE) |>
unite("burrow_species_stage", 'burrow', 'species', 'stage', remove = FALSE) |>
unite("burrow_species", 'burrow', 'species', remove = FALSE)
###1.4.2 Removing redundant BORIS entries
#the next few lines of code are to do with the way BORIS outputs data, where there was redundant empty sections due to both species being independently scored in the software. So here we are filtering out the species specific data from the correlating BORIS species observation and then rebinding the data to remove all of the empty artifact entries.
#extracting just the goby data so that a species variable can be made
goby_dataframe <- decoded_dataframe |>
dplyr::filter(species == "Goby") |>
dplyr::filter(str_detect(video_id, 'g'))
#extracting just the shrimp data so that a species variable can be made
shrimp_dataframe <- decoded_dataframe |>
dplyr::filter(species == "Shrimp") |>
dplyr::filter(str_detect(video_id, 'j'))
#binding the two species specific data frames back together so that there is now no empty
workingdf1 <- bind_rows(goby_dataframe, shrimp_dataframe)
#filtering for the behviours of interest to this experiment
workingdf2 <- workingdf1 |>
group_by(unique_id) |>
dplyr::filter(behaviour == 'INSIDE' | behaviour == 'EDGE' | behaviour == 'NEAR' | behaviour == 'AWAY' |behaviour == 'CONTACT' |behaviour == 'APART' ) |>
mutate(averageduration = sum(duration))
#becuase of the way BORIS outputs data there are redundant individuals in the data, e.g., even if there are only two shrimp in a burrow then BORIS will still output observations for a third shrimp and all value will be 0's. this is what this section of code is removing
workingdf3 <- workingdf2 |>
group_by(burrow, species, individual) |>
dplyr::filter(sum(averageduration) != 0)
###1.4.3 Assigning values to missing entries
#now here we are making the assumption that when a goby or shrimp is present in other videos but missing from the video of focus it is either away and out of frame for the entire trial if it is a goby. Otherwise it is inside the burrow the entire time if it is a shrimp. These assumptions are based off of preliminary observation and the general tendencies of each species. Because of this assumption we are assigning the respective value to the behviour in cases where one of the species is known to be present in a burrow but not viewed for the length of the trial.
fixeddf1 <- workingdf3 |>
dplyr::filter(averageduration == 0) |>
mutate(duration = ifelse(species == "Goby" & behaviour == "AWAY", 300, duration),
duration = ifelse(species == "Shrimp" & behaviour == "INSIDE", 300, duration)
)
fixeddf2 <- workingdf3 |>
dplyr::filter(averageduration != 0)
workingdf4 <- rbind(fixeddf2, fixeddf1) |>
group_by(burrow_species) |>
dplyr::mutate(number_conspecifics = max(as.integer(individual))) |>
ungroup()
#now using the number of conspecifics to calculate a duration value that represents the mean time spent doing that activity for that species (i.e. goby or shrimp) in that burrow as we are un able to preserve identity and thus instead must look at species averages for each behaviour.
workingdf5 <- workingdf4 |>
group_by(burrow_species, behaviour, stage) |>
dplyr::mutate(mean_duration = mean(duration)) |>
ungroup() |>
dplyr::filter(individual == 1)
#transforming the data frame into wide format and turning time spent in each zone into a proportion of the total time (i.e. a value between 0 and 1)
wide_data <- workingdf5 |>
pivot_wider(
id_cols = c(burrow_species, burrow_species_stage, burrow, stage, species, boat_id, stroke, site, number_conspecifics, shrimp_species),
names_from = behaviour,
values_from = mean_duration
)
wide_data <- wide_data |>
mutate(INSIDE = as.numeric(INSIDE),
EDGE = as.numeric(EDGE),
NEAR = as.numeric(NEAR),
AWAY = as.numeric(AWAY),
CONTACT = as.numeric(CONTACT),
APART = as.numeric(APART)) |>
mutate(total = INSIDE+NEAR+EDGE+AWAY,
inside = (INSIDE + EDGE)/(INSIDE+EDGE+NEAR+AWAY),
away = (NEAR + AWAY)/(INSIDE+EDGE+NEAR+AWAY),
number_conspecifics = as.factor(number_conspecifics),
shrimp_contact = CONTACT/(CONTACT+APART)
) |>
dplyr::select(-INSIDE, -EDGE, -NEAR, -AWAY, -CONTACT, -APART)
#checking if every burrow contains observations for both 'Shrimp' and 'Goby'
burrow_values <- unique(wide_data$burrow)
has_shrimp <- burrow_values %in% wide_data$burrow[wide_data$species == 'Shrimp']
has_goby <- burrow_values %in% wide_data$burrow[wide_data$species == 'Goby']
burrows_without_both_species <- burrow_values[!(has_shrimp & has_goby)]
cat("burrows without both 'Shrimp' and 'Goby' species:", paste(burrows_without_both_species, collapse = ', '))
## burrows without both 'Shrimp' and 'Goby' species: 46, 126, 145, 59, 99, 22, 24, 32
#removing the rows with burrows without both 'Shrimp' and 'Goby' species from the dataframe
wide_data <- wide_data[!wide_data$burrow %in% burrows_without_both_species, ]
wide_data <- wide_data |>
mutate(burrow = as.factor(burrow))
#Checking to see the count of observations per burrow
burrow_counts <- table(wide_data$burrow)
print(burrow_counts)
##
## 1 2 3 4 5 6 7 8 9 10 13 15 16 17 18 19 20 21 23 25
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 6
## 26 27 28 29 30 31 33 34 35 36 37 38 39 40 42 44 47 48 49 51
## 6 6 6 6 6 6 5 6 6 6 6 6 6 6 6 6 6 6 6 6
## 52 53 54 55 56 57 58 60 61 62 63 64 65 66 67 68 69 70 71 72
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 93 94 95 96 97 98 100 101 102 103 104 105 106 107 109 110 111 113 114 115
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 116 117 118 119 121 122 123 124 125 127 129 130 131 132 133 134 135 136 137 138
## 6 6 6 5 6 6 6 6 6 6 6 6 6 6 4 6 6 6 6 6
## 139 140 141 142 143 144
## 6 6 5 6 6 6
#there are some burrows that are missing stages due to the animal not being observered at all during that observation so we will make some conditional assumptions and fill in the missing values.
#a function to fill the missing stages for each species within each burrow
fill_missing_stages <- function(df) {
all_stages <- c('Pre', 'During', 'Post')
missing_stages <- setdiff(all_stages, unique(df$stage))
if(length(missing_stages) == 0) return(data.frame())
common_values <- df[1, ]
new_rows <- data.frame()
for(stage in missing_stages) {
new_row <- common_values
new_row$stage <- stage
# Set column values based on species
if(new_row$species == "Goby") {
new_row$away <- 1
new_row$inside <- 0
new_row$total <- 300
new_row$shrimp_contact <- 0
} else if(new_row$species == "Shrimp") {
new_row$away <- 0
new_row$inside <- 1
new_row$total <- 300
new_row$shrimp_contact <- NA
}
new_rows <- rbind(new_rows, new_row)
}
return(new_rows)
}
#applying the function
fixeddf3 <- wide_data |>
group_by(burrow_species) |>
do(fill_missing_stages(.)) |>
ungroup()
wide_data_filled <- rbind(wide_data, fixeddf3)
#Checking to see the count of observations per burrow
filled_burrow_counts <- table(wide_data_filled$burrow)
print(filled_burrow_counts)
##
## 1 2 3 4 5 6 7 8 9 10 13 15 16 17 18 19 20 21 23 25
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 26 27 28 29 30 31 33 34 35 36 37 38 39 40 42 44 47 48 49 51
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 52 53 54 55 56 57 58 60 61 62 63 64 65 66 67 68 69 70 71 72
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 93 94 95 96 97 98 100 101 102 103 104 105 106 107 109 110 111 113 114 115
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 116 117 118 119 121 122 123 124 125 127 129 130 131 132 133 134 135 136 137 138
## 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## 139 140 141 142 143 144
## 6 6 6 6 6 6
###1.4.4 Changing dataframe format
workingdf7 <- wide_data_filled |>
mutate(burrow_stage = paste(burrow, stage))
goby_df <- workingdf7 |>
dplyr::filter(species == "Goby") |>
rename(goby_out = away, goby_total = total, goby_number = number_conspecifics) |>
select(-c("inside", "species", "burrow_species", "burrow_species_stage", "shrimp_contact"))
shrimp_df <- workingdf7 |>
dplyr::filter(species == "Shrimp") |>
rename(shrimp_out = away, shrimp_total = total, shrimp_number = number_conspecifics) |>
select(c("shrimp_out", "shrimp_total", "shrimp_number", "burrow_stage", "shrimp_contact"))
final_df <- inner_join(goby_df, shrimp_df, by = "burrow_stage") |>
mutate(shrimp_contact = ifelse(is.na(shrimp_contact), ifelse(goby_out == 1, 0, 1 - goby_out), shrimp_contact)) |>
filter(site != "SITE O")
write_csv(final_df, "./2-data/final_df.csv")
#2. Statistical analysis using zero-one inflated beta regression ##2.1 Building a global model
#building a global model for the proportion of time gobies spent out of the burrow
zoib_formula <- bf(
# mu (mean) part
goby_out ~ stroke * stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# phi (precision) part
phi ~ stroke * stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# alpha (zero-one-inflation) part
zoi ~ stroke * stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# The one-inflated part, conditional on the 0s
coi ~ stroke * stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
family = zero_one_inflated_beta()
)
##2.2 Specifying priors ###2.2.1 Checking default priors
get_prior(
zoib_formula,
data = final_df
)
## prior class coef group resp dpar
## (flat) b
## (flat) b goby_number2
## (flat) b shrimp_number2
## (flat) b shrimp_number3
## (flat) b shrimp_speciesmannarensis
## (flat) b shrimp_speciessciolii
## (flat) b shrimp_speciesunknown
## (flat) b stagePost
## (flat) b stagePre
## (flat) b stroke4
## (flat) b stroke4:stagePost
## (flat) b stroke4:stagePre
## (flat) b strokecontrol
## (flat) b strokecontrol:stagePost
## (flat) b strokecontrol:stagePre
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd site
## student_t(3, 0, 2.5) sd Intercept site
## student_t(3, 0, 2.5) sd site:burrow
## student_t(3, 0, 2.5) sd Intercept site:burrow
## (flat) b coi
## (flat) b goby_number2 coi
## (flat) b shrimp_number2 coi
## (flat) b shrimp_number3 coi
## (flat) b shrimp_speciesmannarensis coi
## (flat) b shrimp_speciessciolii coi
## (flat) b shrimp_speciesunknown coi
## (flat) b stagePost coi
## (flat) b stagePre coi
## (flat) b stroke4 coi
## (flat) b stroke4:stagePost coi
## (flat) b stroke4:stagePre coi
## (flat) b strokecontrol coi
## (flat) b strokecontrol:stagePost coi
## (flat) b strokecontrol:stagePre coi
## logistic(0, 1) Intercept coi
## student_t(3, 0, 2.5) sd coi
## student_t(3, 0, 2.5) sd site coi
## student_t(3, 0, 2.5) sd Intercept site coi
## student_t(3, 0, 2.5) sd site:burrow coi
## student_t(3, 0, 2.5) sd Intercept site:burrow coi
## (flat) b phi
## (flat) b goby_number2 phi
## (flat) b shrimp_number2 phi
## (flat) b shrimp_number3 phi
## (flat) b shrimp_speciesmannarensis phi
## (flat) b shrimp_speciessciolii phi
## (flat) b shrimp_speciesunknown phi
## (flat) b stagePost phi
## (flat) b stagePre phi
## (flat) b stroke4 phi
## (flat) b stroke4:stagePost phi
## (flat) b stroke4:stagePre phi
## (flat) b strokecontrol phi
## (flat) b strokecontrol:stagePost phi
## (flat) b strokecontrol:stagePre phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi
## student_t(3, 0, 2.5) sd site phi
## student_t(3, 0, 2.5) sd Intercept site phi
## student_t(3, 0, 2.5) sd site:burrow phi
## student_t(3, 0, 2.5) sd Intercept site:burrow phi
## (flat) b zoi
## (flat) b goby_number2 zoi
## (flat) b shrimp_number2 zoi
## (flat) b shrimp_number3 zoi
## (flat) b shrimp_speciesmannarensis zoi
## (flat) b shrimp_speciessciolii zoi
## (flat) b shrimp_speciesunknown zoi
## (flat) b stagePost zoi
## (flat) b stagePre zoi
## (flat) b stroke4 zoi
## (flat) b stroke4:stagePost zoi
## (flat) b stroke4:stagePre zoi
## (flat) b strokecontrol zoi
## (flat) b strokecontrol:stagePost zoi
## (flat) b strokecontrol:stagePre zoi
## logistic(0, 1) Intercept zoi
## student_t(3, 0, 2.5) sd zoi
## student_t(3, 0, 2.5) sd site zoi
## student_t(3, 0, 2.5) sd Intercept site zoi
## student_t(3, 0, 2.5) sd site:burrow zoi
## student_t(3, 0, 2.5) sd Intercept site:burrow zoi
## nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
#we can see that all the beta coefficients are being given a flat prior which means it can take on a value from - infinity to + infinity. All the beta coefficients are on the logit scale so lets plot how changing the value of beta relates to different proportions
###2.2.2 Graphing logit~proportion
tibble(x = seq(-8, 8, by = 0.1)) |>
mutate(y = plogis(x)) |>
ggplot(aes(x = x, y = y)) +
geom_line(size = 1) +
labs(x = "Logit scale", y = "Probability scale")
#we can see that after -5 or +5 the probability doesn't change much, instead of leaving our priors at the default of + and - infinity, we should say that they are most likely to occur between -5 and +5.
###2.2.3 Ploting priors
# Generate data points
x <- seq(-20, 20, length.out = 1000) # Adjusted range for better visualization
# Define parameters for the distributions
loc_logistic <- 0
scale_logistic <- 3
df_student_t <- 3
loc_student_t <- 0
scale_student_t <- 2.5
loc_normal <- 0
scale_normal <- 3
# Calculate PDFs
pdf_logistic <- dlogis(x, location = loc_logistic, scale = scale_logistic)
pdf_student_t <- dt(x, df = df_student_t)
pdf_normal <- dnorm(x, mean = loc_normal, sd = scale_normal)
# Plot the PDFs
df <- data.frame(x = x, Logistic = pdf_logistic, Student_t = pdf_student_t, Normal = pdf_normal)
ggplot(df, aes(x)) +
geom_line(aes(y = Logistic, color = "Logistic(0, 3)")) + # Default priors of Logistic(0,1) have been broadened as these priors seemed too restrictive
geom_line(aes(y = Student_t, color = "Student_t(3, 0, 2.5)")) +
geom_line(aes(y = Normal, color = "Normal(0, 3)")) + #For the beta coefficient prior, we will set it to be normally distributed around 0 with an sd of 5, here is the plotted prior distribution
labs(title = "Probability Density Functions",
x = "x",
y = "Probability Density") +
scale_color_manual(name = "Distribution", values = c("Logistic(0, 3)" = "blue", "Student_t(3, 0, 2.5)" = "green", "Normal(0, 3)" = "red")) +
theme_minimal()
###2.2.4 Setting priors
#the priors for the intercept have been left as the default and the coefficient priors defined above
priors <- c(set_prior("student_t(3, 0, 2.5)", class = "Intercept"),
set_prior("logistic(0, 3)", class = "Intercept", dpar = "coi"),
set_prior("logistic(0, 3)", class = "Intercept", dpar = "zoi"),
set_prior("student_t(3, 0, 2.5)", class = "Intercept", dpar = "phi", lb = 0),
set_prior("student_t(3, 0, 2.5)", class = "sd", lb = 0),
set_prior("student_t(3, 0, 2.5)", class = "sd", dpar = "phi", lb = 0),
set_prior("student_t(3, 0, 2.5)", class = "sd", dpar = "zoi", lb = 0),
set_prior("student_t(3, 0, 2.5)", class = "sd", dpar = "coi", lb = 0),
set_prior("normal(0, 3)", class = "b"),
set_prior("normal(0, 3)", class = "b", dpar = "phi"),
set_prior("normal(0, 3)", class = "b", dpar = "zoi"),
set_prior("normal(0, 3)", class = "b", dpar = "coi"))
###2.2.5 Prior predictive check
prior_formula <- bf(
# mu (mean) part
goby_out ~ stroke*stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# phi (precision) part
phi ~ stroke*stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# alpha (zero-one-inflation) part
zoi ~ stroke*stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# The one-inflated part, conditional on the 0s
coi ~ stroke*stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
family = zero_one_inflated_beta()
)
prior_model <- brm(
prior_formula,
data = final_df,
family = zero_one_inflated_beta(),
prior = priors,
control = list(adapt_delta = 0.99,
max_treedepth = 12),
chains = 4, iter = 1000, warmup = 100,
cores = 4, threads = threading(2),
backend = "cmdstanr",
seed = 1,
sample_prior = "only",
file = "prior_check"
)
pp_check(prior_model, ndraws=20)
pp_check(prior_model, ndraws=100, type = 'stat', stat = 'mean')
##2.3 Modelling the proportion of time gobies spent outside of the burrow
goby_out_formula_1 <- bf(
# mu (mean) part
goby_out ~ stroke*stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# phi (precision) part
phi ~ stroke*stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# alpha (zero-one-inflation) part
zoi ~ stroke*stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
# The one-inflated part, conditional on the 0s
coi ~ stroke*stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow),
family = zero_one_inflated_beta()
)
goby_out_model <- brm(
goby_out_formula_1,
data = final_df,
family = zero_one_inflated_beta(),
prior = priors,
control = list(adapt_delta = 0.99,
max_treedepth = 12),
chains = 4, iter = 4000, warmup = 1000,
cores = 4, threads = threading(2),
backend = "cmdstanr",
seed = 1,
file = "goby_out_model"
)
summary(goby_out_model)
## Family: zero_one_inflated_beta
## Links: mu = logit; phi = log; zoi = logit; coi = logit
## Formula: goby_out ~ stroke * stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow)
## phi ~ stroke * stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow)
## zoi ~ stroke * stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow)
## coi ~ stroke * stage + goby_number + shrimp_number + shrimp_species + (1 | site/burrow)
## Data: final_df (Number of observations: 369)
## Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~site (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.20 0.21 0.01 0.76 1.00 3055 5294
## sd(phi_Intercept) 0.28 0.26 0.01 0.96 1.00 2782 4846
## sd(zoi_Intercept) 0.84 0.60 0.06 2.39 1.00 2557 2582
## sd(coi_Intercept) 1.82 1.70 0.07 6.09 1.00 6696 6013
##
## ~site:burrow (Number of levels: 123)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.95 0.12 0.72 1.21 1.00 2977 5088
## sd(phi_Intercept) 0.18 0.13 0.01 0.49 1.00 2547 4278
## sd(zoi_Intercept) 1.75 0.38 1.07 2.55 1.00 3256 5519
## sd(coi_Intercept) 14.94 8.83 5.25 37.23 1.00 4059 6542
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -0.56 0.58 -1.70 0.58 1.00
## phi_Intercept 0.54 0.56 -0.60 1.61 1.00
## zoi_Intercept -4.65 1.52 -7.77 -1.80 1.00
## coi_Intercept 3.63 4.52 -4.69 13.29 1.00
## stroke4 -0.43 0.34 -1.09 0.23 1.00
## strokecontrol 0.25 0.33 -0.41 0.90 1.00
## stagePost -0.04 0.23 -0.49 0.40 1.00
## stagePre 0.31 0.23 -0.14 0.76 1.00
## goby_number2 0.76 0.27 0.23 1.29 1.00
## shrimp_number2 0.69 0.28 0.14 1.25 1.00
## shrimp_number3 0.38 0.84 -1.30 2.05 1.00
## shrimp_speciesmannarensis 0.41 0.50 -0.58 1.37 1.00
## shrimp_speciessciolii 0.31 0.69 -1.06 1.68 1.00
## shrimp_speciesunknown 0.47 0.55 -0.60 1.54 1.00
## stroke4:stagePost 0.37 0.36 -0.36 1.07 1.00
## strokecontrol:stagePost 0.11 0.38 -0.64 0.86 1.00
## stroke4:stagePre 0.49 0.37 -0.25 1.22 1.00
## strokecontrol:stagePre -0.60 0.36 -1.31 0.11 1.00
## phi_stroke4 -0.75 0.38 -1.50 -0.01 1.00
## phi_strokecontrol -0.13 0.43 -0.97 0.71 1.00
## phi_stagePost -0.40 0.36 -1.13 0.30 1.00
## phi_stagePre -0.22 0.37 -0.96 0.48 1.00
## phi_goby_number2 0.97 0.27 0.45 1.49 1.00
## phi_shrimp_number2 0.22 0.25 -0.26 0.73 1.00
## phi_shrimp_number3 -0.16 0.72 -1.55 1.25 1.00
## phi_shrimp_speciesmannarensis 0.58 0.42 -0.21 1.46 1.00
## phi_shrimp_speciessciolii 0.20 0.55 -0.87 1.30 1.00
## phi_shrimp_speciesunknown 0.49 0.48 -0.42 1.46 1.00
## phi_stroke4:stagePost 0.63 0.51 -0.35 1.62 1.00
## phi_strokecontrol:stagePost 0.08 0.56 -1.00 1.16 1.00
## phi_stroke4:stagePre 0.53 0.54 -0.53 1.57 1.00
## phi_strokecontrol:stagePre -0.06 0.57 -1.17 1.06 1.00
## zoi_stroke4 0.20 0.83 -1.42 1.84 1.00
## zoi_strokecontrol 1.54 0.77 0.08 3.09 1.00
## zoi_stagePost -0.03 0.67 -1.33 1.27 1.00
## zoi_stagePre 0.99 0.63 -0.23 2.25 1.00
## zoi_goby_number2 -1.57 0.71 -3.05 -0.23 1.00
## zoi_shrimp_number2 0.02 0.61 -1.17 1.24 1.00
## zoi_shrimp_number3 0.44 1.50 -2.57 3.33 1.00
## zoi_shrimp_speciesmannarensis 1.89 1.29 -0.55 4.56 1.00
## zoi_shrimp_speciessciolii 1.97 1.47 -0.87 4.90 1.00
## zoi_shrimp_speciesunknown 2.28 1.34 -0.26 4.99 1.00
## zoi_stroke4:stagePost 1.28 0.94 -0.53 3.15 1.00
## zoi_strokecontrol:stagePost 0.45 0.88 -1.27 2.20 1.00
## zoi_stroke4:stagePre -0.41 0.93 -2.24 1.43 1.00
## zoi_strokecontrol:stagePre -1.13 0.87 -2.83 0.55 1.00
## coi_stroke4 -0.97 2.58 -5.87 4.22 1.00
## coi_strokecontrol -0.64 2.62 -5.62 4.77 1.00
## coi_stagePost 0.59 2.02 -3.27 4.74 1.00
## coi_stagePre -1.34 2.01 -5.50 2.44 1.00
## coi_goby_number2 1.87 2.89 -3.88 7.46 1.00
## coi_shrimp_number2 1.88 2.65 -3.43 6.98 1.00
## coi_shrimp_number3 0.58 2.93 -5.23 6.29 1.00
## coi_shrimp_speciesmannarensis -0.18 2.61 -5.22 5.00 1.00
## coi_shrimp_speciessciolii 0.11 2.84 -5.52 5.63 1.00
## coi_shrimp_speciesunknown 0.16 2.67 -5.09 5.36 1.00
## coi_stroke4:stagePost -0.67 2.64 -5.78 4.52 1.00
## coi_strokecontrol:stagePost 0.07 2.52 -4.94 5.10 1.00
## coi_stroke4:stagePre 0.60 2.68 -4.63 5.80 1.00
## coi_strokecontrol:stagePre -0.34 2.55 -5.34 4.75 1.00
## Bulk_ESS Tail_ESS
## Intercept 3817 6057
## phi_Intercept 3954 5442
## zoi_Intercept 4911 6500
## coi_Intercept 7410 7550
## stroke4 4145 6511
## strokecontrol 4019 6210
## stagePost 5396 7566
## stagePre 5218 7587
## goby_number2 3788 6321
## shrimp_number2 3977 6505
## shrimp_number3 6113 7672
## shrimp_speciesmannarensis 3627 5992
## shrimp_speciessciolii 4123 6904
## shrimp_speciesunknown 3833 6301
## stroke4:stagePost 5684 7467
## strokecontrol:stagePost 5607 7906
## stroke4:stagePre 5496 7223
## strokecontrol:stagePre 5464 7712
## phi_stroke4 3952 5841
## phi_strokecontrol 3759 5665
## phi_stagePost 3835 5937
## phi_stagePre 4177 6664
## phi_goby_number2 7361 8798
## phi_shrimp_number2 7372 8376
## phi_shrimp_number3 8680 8719
## phi_shrimp_speciesmannarensis 4476 6512
## phi_shrimp_speciessciolii 5578 7638
## phi_shrimp_speciesunknown 4466 6560
## phi_stroke4:stagePost 4803 6699
## phi_strokecontrol:stagePost 4285 6916
## phi_stroke4:stagePre 4726 6846
## phi_strokecontrol:stagePre 4172 6623
## zoi_stroke4 4669 7059
## zoi_strokecontrol 4447 6447
## zoi_stagePost 4595 7525
## zoi_stagePre 5029 7451
## zoi_goby_number2 6448 7257
## zoi_shrimp_number2 5830 7445
## zoi_shrimp_number3 8548 9025
## zoi_shrimp_speciesmannarensis 5219 6637
## zoi_shrimp_speciessciolii 5589 7427
## zoi_shrimp_speciesunknown 5358 7012
## zoi_stroke4:stagePost 5029 7911
## zoi_strokecontrol:stagePost 5072 7101
## zoi_stroke4:stagePre 5527 8009
## zoi_strokecontrol:stagePre 5585 7375
## coi_stroke4 9982 8939
## coi_strokecontrol 8616 8392
## coi_stagePost 9837 8241
## coi_stagePre 9290 8901
## coi_goby_number2 9040 8784
## coi_shrimp_number2 8762 8986
## coi_shrimp_number3 13937 8873
## coi_shrimp_speciesmannarensis 11322 9704
## coi_shrimp_speciessciolii 14206 9008
## coi_shrimp_speciesunknown 10515 9219
## coi_stroke4:stagePost 11886 9584
## coi_strokecontrol:stagePost 11208 9021
## coi_stroke4:stagePre 12132 8898
## coi_strokecontrol:stagePre 12755 8356
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
###2.3.1 Model fit and diagnostics ####2.3.1.1 Model convergence
#convergence diagnostics (R-hat)
extract_rhats <- function(model) {
model_summary <- summary(model)
#extract R-hat values for fixed effects
fixed_rhats <- model_summary$fixed[, "Rhat", drop = FALSE]
#check for random effects and extract R-hat values if present
if (!is.null(model_summary$random)) {
random_rhats <- do.call(rbind, lapply(model_summary$random, function(x) x[, "Rhat", drop = FALSE]))
#combine fixed and random effects R-hat values
rhats <- rbind(fixed_rhats, random_rhats)
} else {
rhats <- fixed_rhats
}
return(rhats)
}
extract_rhats(goby_out_model)
## Rhat
## Intercept 1.0005578
## phi_Intercept 1.0009595
## zoi_Intercept 0.9999930
## coi_Intercept 1.0002505
## stroke4 1.0002134
## strokecontrol 1.0000376
## stagePost 1.0007459
## stagePre 1.0012678
## goby_number2 1.0002532
## shrimp_number2 1.0010880
## shrimp_number3 1.0006650
## shrimp_speciesmannarensis 1.0013751
## shrimp_speciessciolii 1.0010367
## shrimp_speciesunknown 1.0011741
## stroke4:stagePost 1.0006817
## strokecontrol:stagePost 1.0003990
## stroke4:stagePre 1.0006277
## strokecontrol:stagePre 1.0008589
## phi_stroke4 1.0005548
## phi_strokecontrol 1.0000521
## phi_stagePost 1.0003700
## phi_stagePre 1.0011527
## phi_goby_number2 1.0005467
## phi_shrimp_number2 1.0001211
## phi_shrimp_number3 0.9999848
## phi_shrimp_speciesmannarensis 1.0010171
## phi_shrimp_speciessciolii 1.0007876
## phi_shrimp_speciesunknown 1.0005111
## phi_stroke4:stagePost 1.0005792
## phi_strokecontrol:stagePost 0.9999932
## phi_stroke4:stagePre 1.0002913
## phi_strokecontrol:stagePre 1.0001436
## zoi_stroke4 1.0003714
## zoi_strokecontrol 1.0010740
## zoi_stagePost 1.0000282
## zoi_stagePre 1.0004724
## zoi_goby_number2 1.0006268
## zoi_shrimp_number2 1.0000590
## zoi_shrimp_number3 1.0000891
## zoi_shrimp_speciesmannarensis 1.0005980
## zoi_shrimp_speciessciolii 1.0006451
## zoi_shrimp_speciesunknown 1.0004113
## zoi_stroke4:stagePost 0.9999697
## zoi_strokecontrol:stagePost 1.0002877
## zoi_stroke4:stagePre 1.0006642
## zoi_strokecontrol:stagePre 1.0003515
## coi_stroke4 1.0001140
## coi_strokecontrol 1.0001493
## coi_stagePost 1.0001532
## coi_stagePre 1.0002669
## coi_goby_number2 1.0007911
## coi_shrimp_number2 1.0003449
## coi_shrimp_number3 1.0000471
## coi_shrimp_speciesmannarensis 1.0001202
## coi_shrimp_speciessciolii 0.9999673
## coi_shrimp_speciesunknown 1.0002196
## coi_stroke4:stagePost 1.0000630
## coi_strokecontrol:stagePost 1.0000740
## coi_stroke4:stagePre 1.0005086
## coi_strokecontrol:stagePre 1.0002274
## site.sd(Intercept) 1.0011060
## site.sd(phi_Intercept) 1.0006338
## site.sd(zoi_Intercept) 1.0006809
## site.sd(coi_Intercept) 1.0006903
## site:burrow.sd(Intercept) 1.0000363
## site:burrow.sd(phi_Intercept) 1.0007214
## site:burrow.sd(zoi_Intercept) 1.0010080
## site:burrow.sd(coi_Intercept) 1.0009396
####2.3.1.2 Model fit
#posterior predictive checks
pp_check(goby_out_model, ndraws = 50)
pp_check(goby_out_model, type='stat', stat='mean')
#fitted vs predicted values
summary(fitted(goby_out_model))
## Estimate Est.Error Q2.5 Q97.5
## Min. :0.1588 Min. :0.04909 Min. :0.007135 Min. :0.4130
## 1st Qu.:0.4984 1st Qu.:0.10560 1st Qu.:0.213903 1st Qu.:0.7719
## Median :0.6774 Median :0.12828 Median :0.379441 Median :0.8931
## Mean :0.6297 Mean :0.12544 Mean :0.367441 Mean :0.8473
## 3rd Qu.:0.7918 3rd Qu.:0.14591 3rd Qu.:0.516355 3rd Qu.:0.9451
## Max. :0.9326 Max. :0.24505 Max. :0.803405 Max. :0.9987
summary(predict(goby_out_model))
## Estimate Est.Error Q2.5 Q97.5
## Min. :0.1595 Min. :0.1284 Min. :0.00000 Min. :0.9330
## 1st Qu.:0.5003 1st Qu.:0.2471 1st Qu.:0.00000 1st Qu.:1.0000
## Median :0.6779 Median :0.2965 Median :0.00000 Median :1.0000
## Mean :0.6294 Mean :0.2854 Mean :0.06544 Mean :0.9992
## 3rd Qu.:0.7922 3rd Qu.:0.3295 3rd Qu.:0.04008 3rd Qu.:1.0000
## Max. :0.9328 Max. :0.4557 Max. :0.57493 Max. :1.0000
###2.3.2 Posterior marginal effects ####2.3.2.1 Pairwise contrasts of stroke and stage
goby_out_contrast_effect <- goby_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
# Extract the relevant contrast information
goby_out_results.table <- as.data.frame(goby_out_contrast_effect) |>
dplyr::filter(contrast %in% c("control Pre - 2 Pre", "control Pre - 4 Pre","4 Pre - 2 Pre", "control Pre - control During", "control Pre - control Post", "4 Pre - 4 During", "4 Pre - 4 Post", "2 Pre - 2 During", "2 Pre - 2 Post"))
goby_out_results.table
## contrast estimate lower.HPD upper.HPD
## 2 Pre - 2 During 0.06131526 -0.04085772 0.16578566
## 2 Pre - 2 Post 0.06697584 -0.04508346 0.18552813
## 4 Pre - 4 During 0.15184351 0.01311128 0.28627348
## 4 Pre - 4 Post 0.06674910 -0.08135213 0.22150063
## 4 Pre - 2 Pre 0.00457735 -0.15864033 0.13777826
## control Pre - control During -0.06948150 -0.20576921 0.06260540
## control Pre - control Post -0.09180295 -0.25366372 0.05686286
## control Pre - 2 Pre -0.06730160 -0.22705272 0.08824646
## control Pre - 4 Pre -0.07129914 -0.25472546 0.10683535
##
## Results are averaged over the levels of: goby_number, shrimp_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
goby_out_contrast_draws <- goby_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise") |>
gather_emmeans_draws()
goby_out_contrast_draws
## # A tibble: 432,000 × 5
## # Groups: contrast [36]
## contrast .chain .iteration .draw .value
## <chr> <int> <int> <int> <dbl>
## 1 4 During - 2 During NA NA 1 -0.00820
## 2 4 During - 2 During NA NA 2 -0.0491
## 3 4 During - 2 During NA NA 3 -0.0492
## 4 4 During - 2 During NA NA 4 -0.107
## 5 4 During - 2 During NA NA 5 -0.0321
## 6 4 During - 2 During NA NA 6 -0.0724
## 7 4 During - 2 During NA NA 7 -0.138
## 8 4 During - 2 During NA NA 8 -0.0889
## 9 4 During - 2 During NA NA 9 -0.0438
## 10 4 During - 2 During NA NA 10 -0.0502
## # ℹ 431,990 more rows
goby_out_temporal_contrasts <- goby_out_contrast_draws |>
dplyr::filter(contrast %in% c("control Pre - control During", "control Pre - control Post", "4 Pre - 4 During", "4 Pre - 4 Post", "2 Pre - 2 During", "2 Pre - 2 Post")) |>
mutate(stroke = ifelse(contrast %in% c("control Pre - control During", "control Pre - control Post"), "Control", ifelse(contrast %in% c("4 Pre - 4 During", "4 Pre - 4 Post"), "4", "2")))
goby_out_temporal_contrasts$stroke <- factor(goby_out_temporal_contrasts$stroke, levels = c("Control", "4", "2"))
fig.goby_out_temporal_contrasts<- ggplot(goby_out_temporal_contrasts, aes(x = .value, y = factor(contrast, levels = c("control Pre - control During", "control Pre - control Post", "4 Pre - 4 During", "4 Pre - 4 Post", "2 Pre - 2 During", "2 Pre - 2 Post")), color = stroke, fill = stroke)) +
geom_vline(xintercept = 0, color = "black", size = 0.5) +
stat_halfeye(.width = c(0.89, 0.95))+
scale_color_manual(values = c(colour11, colour12, colour13)) +
scale_fill_manual(values = alpha(c(colour1, colour2, colour3), 0.8)) +
labs(y = "Contrasts of the time periods relative to noise exposure",
fill = "Noise treatment", color = "Noise treatment",
fill_ramp = "Credible interval") +
xlab(NULL) +
coord_flip() +
theme(legend.position = "right",
axis.text.x = element_text(angle = 50, hjust =1)) +
scale_y_discrete(labels = c("Pre - During", "Pre - Post", "Pre - During", "Pre - Post", "Pre - During", "Pre - Post")) +
scale_x_continuous(breaks = seq(-0.4, 0.5, 0.1), labels = c("-0.4","-0.3", "-0.2","-0.1" , "0", "0.1", "0.2", "0.3", "0.4", "0.5"), limits = c(-0.4, 0.4))
fig.goby_out_temporal_contrasts
###2.3.3 Extracting and plotting model estimates ####2.3.3.1 Extracting
estimates
goby_out_model_estimates <- goby_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA)
goby_out_model_estimates <- as.data.frame(goby_out_model_estimates)
#Storing the model estimate means for each group to plot later on
goby_out_model_estimates_control <- goby_out_model_estimates |>
dplyr::filter(stroke == "control") |>
rename(.value = emmean)
goby_out_model_estimates_four <- goby_out_model_estimates |>
dplyr::filter(stroke == "4")|>
rename(.value = emmean)
goby_out_model_estimates_two <- goby_out_model_estimates |>
dplyr::filter(stroke == "2")|>
rename(.value = emmean)
####2.3.3.2 Plotting model estimates and raw data
goby_out_model_draws <- goby_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA) |>
gather_emmeans_draws()
goby_out_model_draws <- goby_out_model_draws |>
rename(emmean = .value)
goby_out_raw_data <- final_df |>
rename(emmean = goby_out)
goby_out_raw_data$stage <- factor(goby_out_raw_data$stage, levels = c("Pre", "During", "Post"))
goby_out_raw_data$stroke <- factor(goby_out_raw_data$stroke, levels = c("control", "4", "2"))
goby_out_model_estimates$stage <- factor(goby_out_model_estimates$stage, levels = c("Pre", "During", "Post"))
goby_out_model_estimates$stroke <- factor(goby_out_model_estimates$stroke, levels = c("control", "4", "2"))
goby_out_model_draws$stage <- factor(goby_out_model_draws$stage, levels = c("Pre", "During", "Post"))
goby_out_model_draws$stroke <- factor(goby_out_model_draws$stroke, levels = c("control", "4", "2"))
goby_out_means <- goby_out_raw_data |>
group_by(stage, stroke) |>
summarise(emmean = mean(emmean))
goby_out_fig_together <- ggplot(data = goby_out_raw_data, aes(x = stage, y = emmean, color = stroke, fill = stroke)) +
geom_point(data = goby_out_means, alpha = 1, position = position_dodgenudge(width = 0.8, x = -0.09), size = 9, shape = "_") +
geom_half_point(alpha = 0.3, position = position_dodgenudge(width = 0.8, x = -0.03), side = "l", range_scale = 0.5) +
geom_half_violin(data = goby_out_model_draws, aes(x = stage, y = emmean, color = stroke, fill = stroke), alpha = 0.6, position = position_dodge(width = 0.8), side = "r", size = 0) +
geom_point(data = goby_out_model_estimates, aes(x = stage, y = emmean, color = stroke), position = position_dodge(width = 0.8), size = 3) +
geom_errorbar(data = goby_out_model_estimates, aes(x = stage, y = emmean, ymin = lower.HPD, ymax = upper.HPD, color = stroke), position = position_dodge(width = 0.8), width = 0, size = 1.25) +
scale_fill_manual(values = c("control" = "#76A34A", "4" = "#F29E00", "2" = "#E25A00"), name = "Noise exposure", guide = guide_legend(override.aes = list(size = 3))) +
scale_color_manual(values = c("control" = "#2A330E", "4" = "#805213", "2" = "#542709"), name = "Noise exposure", guide = guide_legend(override.aes = list(size = 3))) +
scale_x_discrete(labels = c("Pre", "During", "Post")) +
labs(title = "The effects of boat noise on goby burrow use",
y = "Proportion of time in out of the burrow",
x = "Experimental phase relative to noise exposure") +
theme_classic() +
theme(
axis.text.x = element_text(size = 10, margin = margin(t = 5)),
axis.text.y = element_text(size = 10, margin = margin(r = 5)),
axis.title.x = element_text(size = 13, margin = margin(t = 10)),
axis.title.y = element_text(size = 13, margin = margin(r = 10)),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
legend.position = "right",
legend.justification = "top",
plot.title = element_text(hjust = 0.5, size = 16),
plot.margin = margin(10, 10, 10, 10)
)
goby_out_fig_together
####2.3.3.3 Plotting control estimates and raw data
goby_out_model_draws <- goby_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA) |>
gather_emmeans_draws()
goby_out_raw_data <- final_df |>
rename(.value = goby_out)
goby_out_model_draws$stage <- factor(goby_out_model_draws$stage, levels = c("Pre", "During", "Post"))
goby_out_model_draws$stroke <- factor(goby_out_model_draws$stroke, levels = c("control", "4", "2"))
goby_out_control_model_draws <- goby_out_model_draws |>
dplyr::filter(stroke == "control")
goby_control_raw_data <- final_df |>
dplyr::filter(stroke =="control") |>
rename(.value = goby_out)
goby_control_means <- goby_control_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
goby_control_raw_data$stage <- factor(goby_control_raw_data$stage, levels = c("Pre", "During", "Post"))
#assigning a value for the baseline mean
goby_out_model_estimates_control_intercept <- goby_out_model_estimates_control %>%
dplyr::filter(stage == "Pre") %>%
pull(.value)
goby_out_fig_control <- ggplot(goby_control_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = goby_out_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = goby_out_model_estimates_control, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = goby_out_model_estimates_control, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = goby_control_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
labs(x="Period relative to treatment",
y="Proportion of time out of burrow") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_text(size = 13),
axis.text.y = element_text(size = 13),
panel.grid = element_blank(),
axis.title.y = element_text(size = 16, face = "bold", margin = margin(t = 0, b = 0, l = 0, r = 10)),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#A3C088", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("Control")+
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#76A34A", "Post" = "#D3D9A7")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#2A330E", "Post" = "#5d634c"))
goby_out_fig_control
####2.3.3.4 Plotting 4-stroke estimates and raw data
goby_out_control_model_draws <- goby_out_model_draws |>
dplyr::filter(stroke == "4")
goby_four_raw_data <- final_df |>
dplyr::filter(stroke =="4") |>
rename(.value = goby_out)
goby_four_raw_data$stage <- factor(goby_four_raw_data$stage, levels = c("Pre", "During", "Post"))
goby_four_means <- goby_four_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
goby_out_fig_four <- ggplot(goby_four_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = goby_out_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = goby_out_model_estimates_four, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = goby_out_model_estimates_four, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = goby_four_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
xlab("Experimental phase relative to noise exposure")+
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.x = element_text(size = 13),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid = element_blank(),
axis.title.x = element_text(size = 16, face = "bold", margin = margin(t = 15, b = 0, l = 0, r = 0)),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#F6BC65", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("4-stroke")+
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#F29E00", "Post" = "#F2D5A0")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#805213", "Post" = "#70532b"))
goby_out_fig_four
####2.3.3.5 Plotting 2-stroke estimates and raw data
goby_out_control_model_draws <- goby_out_model_draws |>
dplyr::filter(stroke == "2")
goby_two_raw_data <- final_df |>
dplyr::filter(stroke =="2") |>
rename(.value = goby_out)
goby_two_raw_data$stage <- factor(goby_two_raw_data$stage, levels = c("Pre", "During", "Post"))
goby_two_means <- goby_two_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
goby_out_fig_two <- ggplot(goby_two_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = goby_out_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = goby_out_model_estimates_two, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = goby_out_model_estimates_two, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = goby_two_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_text(size = 13),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid = element_blank(),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#EB9063", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("2-stroke") +
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#E25A00", "Post" = "#F2B389")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#542709", "Post" = "#66422a"))
goby_out_fig_two
####2.3.3.6 Combining plots into one figure
goby_out_plot3 <- goby_out_fig_control + goby_out_fig_four + goby_out_fig_two +
plot_layout(ncol = 3)
goby_out_plot3
ggsave("./3-figs/goby_out_plot3.jpeg", plot = goby_out_plot3, units = "px", width = 2200, height = 1400, dpi = 300)
###2.3.4 Effects of covariates on the proportion of time gobies spent outside of the burrow ####2.3.4.1 Goby number
goby_out_model |>
emmeans(~ goby_number,
epred = TRUE,
re_formla = NULL)
## goby_number emmean lower.HPD upper.HPD
## 1 0.575 0.423 0.729
## 2 0.707 0.549 0.842
##
## Results are averaged over the levels of: stroke, stage, shrimp_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
goby_out_goby_number <- goby_out_model |>
emmeans(~ goby_number,
epred = TRUE,
re_formla = NULL) |>
contrast(method = "revpairwise")
goby_out_goby_number
## contrast estimate lower.HPD upper.HPD
## goby_number2 - goby_number1 0.128 -0.00135 0.249
##
## Results are averaged over the levels of: stroke, stage, shrimp_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
goby_out_goby_number_draws <- goby_out_model |>
emmeans(~ goby_number,
epred = TRUE,
re_formla = NULL) |>
contrast(method = "revpairwise")|>
gather_emmeans_draws()
goby_out_goby_number_fig<- ggplot(goby_out_goby_number_draws, aes(x = .value, y = contrast)) +
geom_vline(xintercept = 0, color = "black", size = 0.5) +
stat_halfeye(.width = c(0.95))+
labs(y = "Contrasts of the effect of the number of goby",
fill = "Noise treatment", color = "Noise treatment",
fill_ramp = "Credible interval") +
xlab(NULL) +
coord_flip()
goby_out_goby_number_fig
####2.3.4.2 Shrimp number
goby_out_shrimp_number <- goby_out_model |>
emmeans(~ shrimp_number,
epred = TRUE,
re_formla = NULL) |>
contrast(method = "revpairwise")
goby_out_shrimp_number
## contrast estimate lower.HPD upper.HPD
## shrimp_number2 - shrimp_number1 0.1431 0.025 0.269
## shrimp_number3 - shrimp_number1 0.0867 -0.218 0.355
## shrimp_number3 - shrimp_number2 -0.0574 -0.366 0.187
##
## Results are averaged over the levels of: stroke, stage, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
goby_out_shrimp_number_estimates <- goby_out_model |>
emmeans(~ shrimp_number,
epred = TRUE,
re_formla = NULL)
goby_out_shrimp_number_estimates
## shrimp_number emmean lower.HPD upper.HPD
## 1 0.566 0.413 0.708
## 2 0.711 0.614 0.803
## 3 0.651 0.339 0.909
##
## Results are averaged over the levels of: stroke, stage, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
goby_out_shrimp_number_draws <- goby_out_model |>
emmeans(~ shrimp_number,
epred = TRUE,
re_formla = NULL) |>
contrast(method = "revpairwise")|>
gather_emmeans_draws()
goby_out_shrimp_number_fig<- ggplot(goby_out_shrimp_number_draws, aes(x = .value, y = contrast)) +
geom_vline(xintercept = 0, color = "black", size = 0.5) +
stat_halfeye(.width = c(0.95))+
labs(y = "Contrasts of the effect of the number of shrimp",
fill = "Noise treatment", color = "Noise treatment",
fill_ramp = "Credible interval") +
xlab(NULL) +
coord_flip()
goby_out_shrimp_number_fig
####2.3.4.3 Shrimp species
goby_out_shrimp_species <- goby_out_model |>
emmeans(~ shrimp_species,
epred = TRUE,
re_formla = NULL) |>
contrast(method = "revpairwise")
goby_out_shrimp_species <- as.data.frame(goby_out_shrimp_species)
goby_out_shrimp_species_estimates <- goby_out_model |>
emmeans(~ shrimp_species,
epred = TRUE,
re_formla = NULL)
goby_out_shrimp_species_estimates <- as.data.frame(goby_out_shrimp_species_estimates)
goby_out_shrimp_species_draws <- goby_out_model |>
emmeans(~ shrimp_species,
epred = TRUE,
re_formla = NULL) |>
contrast(method = "revpairwise")|>
gather_emmeans_draws()
goby_out_shrimp_species_fig<- ggplot(goby_out_shrimp_species_draws, aes(x = .value, y = contrast)) +
geom_vline(xintercept = 0, color = "black", size = 0.5) +
stat_halfeye(.width = c(0.95))+
labs(y = "Contrasts of the effect of the species of shrimp",
fill = "Noise treatment", color = "Noise treatment",
fill_ramp = "Credible interval") +
xlab(NULL) +
coord_flip()
goby_out_shrimp_species_fig
##2.4 Modelling the proportion of time shrimp spent outside of the burrow
shrimp_out_formula <- bf(
# mu (mean) part
shrimp_out ~ stroke*stage + shrimp_number + goby_number + shrimp_species + (1|site/burrow),
# phi (precision) part
phi ~ stroke*stage + shrimp_number + goby_number + shrimp_species + (1|site/burrow),
# alpha (zero-one-inflation) part
zoi ~ stroke*stage + shrimp_number + goby_number + shrimp_species + (1|site/burrow),
# The one-inflated part, conditional on the 0s
coi ~ stroke*stage + shrimp_number + goby_number + shrimp_species + (1|site/burrow),
family = zero_one_inflated_beta()
)
shrimp_out_model <- brm(
shrimp_out_formula,
data = final_df,
family = zero_one_inflated_beta(),
prior = priors,
control = list(adapt_delta = 0.99,
max_treedepth = 12),
chains = 4, iter = 4000, warmup = 1000,
cores = 4, threads = threading(2),
backend = "cmdstanr",
seed = 1,
file = "shrimp_out_model"
)
summary(shrimp_out_model)
## Family: zero_one_inflated_beta
## Links: mu = logit; phi = log; zoi = logit; coi = logit
## Formula: shrimp_out ~ stroke * stage + shrimp_number + goby_number + shrimp_species + (1 | site/burrow)
## phi ~ stroke * stage + shrimp_number + goby_number + shrimp_species + (1 | site/burrow)
## zoi ~ stroke * stage + shrimp_number + goby_number + shrimp_species + (1 | site/burrow)
## coi ~ stroke * stage + shrimp_number + goby_number + shrimp_species + (1 | site/burrow)
## Data: final_df (Number of observations: 369)
## Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~site (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.18 0.18 0.00 0.67 1.00 3882 5427
## sd(phi_Intercept) 0.28 0.26 0.01 0.97 1.00 2972 5304
## sd(zoi_Intercept) 0.53 0.51 0.02 1.84 1.00 5979 6908
## sd(coi_Intercept) 2.26 2.34 0.08 7.93 1.00 9225 7074
##
## ~site:burrow (Number of levels: 123)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.77 0.10 0.58 0.96 1.00 3318 6077
## sd(phi_Intercept) 0.32 0.18 0.02 0.68 1.00 1797 3263
## sd(zoi_Intercept) 2.14 0.60 1.10 3.46 1.00 3039 4614
## sd(coi_Intercept) 2.04 1.77 0.08 6.46 1.00 7376 7213
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -1.63 0.50 -2.60 -0.62 1.00
## phi_Intercept 1.60 0.63 0.32 2.79 1.00
## zoi_Intercept -2.24 1.49 -5.36 0.60 1.00
## coi_Intercept -7.26 4.01 -15.84 -0.13 1.00
## stroke4 0.21 0.27 -0.30 0.74 1.00
## strokecontrol 0.38 0.27 -0.15 0.91 1.00
## stagePost 0.23 0.19 -0.15 0.61 1.00
## stagePre 0.71 0.19 0.34 1.09 1.00
## shrimp_number2 0.52 0.26 0.01 1.02 1.00
## shrimp_number3 0.46 0.54 -0.63 1.54 1.00
## goby_number2 0.35 0.22 -0.09 0.79 1.00
## shrimp_speciesmannarensis 0.01 0.42 -0.85 0.82 1.00
## shrimp_speciessciolii -0.61 0.54 -1.68 0.43 1.00
## shrimp_speciesunknown -0.04 0.46 -0.96 0.85 1.00
## stroke4:stagePost 0.17 0.29 -0.41 0.74 1.00
## strokecontrol:stagePost 0.00 0.28 -0.55 0.55 1.00
## stroke4:stagePre -0.23 0.30 -0.84 0.36 1.00
## strokecontrol:stagePre -0.57 0.31 -1.19 0.03 1.00
## phi_stroke4 -0.31 0.46 -1.20 0.62 1.00
## phi_strokecontrol -0.61 0.42 -1.43 0.22 1.00
## phi_stagePost -0.55 0.39 -1.31 0.21 1.00
## phi_stagePre -0.58 0.43 -1.42 0.28 1.00
## phi_shrimp_number2 0.74 0.26 0.24 1.24 1.00
## phi_shrimp_number3 2.11 0.78 0.56 3.62 1.00
## phi_goby_number2 0.06 0.25 -0.44 0.55 1.00
## phi_shrimp_speciesmannarensis -0.06 0.50 -1.00 0.94 1.00
## phi_shrimp_speciessciolii 0.15 0.57 -0.98 1.26 1.00
## phi_shrimp_speciesunknown 0.07 0.53 -0.92 1.16 1.00
## phi_stroke4:stagePost -0.15 0.56 -1.25 0.95 1.00
## phi_strokecontrol:stagePost 0.33 0.54 -0.75 1.39 1.00
## phi_stroke4:stagePre -0.24 0.62 -1.47 0.96 1.00
## phi_strokecontrol:stagePre 0.23 0.58 -0.90 1.36 1.00
## zoi_stroke4 -0.99 1.01 -3.01 0.90 1.00
## zoi_strokecontrol -1.31 1.05 -3.48 0.70 1.00
## zoi_stagePost 0.30 0.67 -1.00 1.60 1.00
## zoi_stagePre 0.31 0.67 -1.00 1.64 1.00
## zoi_shrimp_number2 -1.89 0.76 -3.50 -0.54 1.00
## zoi_shrimp_number3 -0.75 1.73 -4.26 2.50 1.00
## zoi_goby_number2 -0.41 0.82 -2.08 1.18 1.00
## zoi_shrimp_speciesmannarensis 0.53 1.30 -1.95 3.16 1.00
## zoi_shrimp_speciessciolii -0.46 1.75 -4.00 2.86 1.00
## zoi_shrimp_speciesunknown 0.52 1.40 -2.24 3.35 1.00
## zoi_stroke4:stagePost 0.16 1.15 -2.08 2.40 1.00
## zoi_strokecontrol:stagePost -0.36 1.24 -2.85 2.09 1.00
## zoi_stroke4:stagePre 0.59 1.10 -1.53 2.82 1.00
## zoi_strokecontrol:stagePre -0.35 1.21 -2.78 1.96 1.00
## coi_stroke4 2.82 2.39 -1.91 7.61 1.00
## coi_strokecontrol -0.35 2.64 -5.65 4.74 1.00
## coi_stagePost -1.42 2.45 -6.30 3.24 1.00
## coi_stagePre -1.73 2.38 -6.54 2.83 1.00
## coi_shrimp_number2 -1.51 2.35 -6.10 3.11 1.00
## coi_shrimp_number3 -0.17 2.92 -5.89 5.42 1.00
## coi_goby_number2 -0.96 2.59 -6.06 4.12 1.00
## coi_shrimp_speciesmannarensis 0.92 2.59 -4.09 6.06 1.00
## coi_shrimp_speciessciolii -0.21 2.88 -5.92 5.41 1.00
## coi_shrimp_speciesunknown -0.76 2.68 -6.11 4.44 1.00
## coi_stroke4:stagePost -0.68 2.74 -6.12 4.63 1.00
## coi_strokecontrol:stagePost 0.02 2.85 -5.58 5.53 1.00
## coi_stroke4:stagePre -0.82 2.70 -6.16 4.36 1.00
## coi_strokecontrol:stagePre -0.06 2.85 -5.68 5.37 1.00
## Bulk_ESS Tail_ESS
## Intercept 6281 7699
## phi_Intercept 5987 6801
## zoi_Intercept 9282 8789
## coi_Intercept 12638 9021
## stroke4 7249 8499
## strokecontrol 6875 8837
## stagePost 7761 9210
## stagePre 7801 8761
## shrimp_number2 8817 9194
## shrimp_number3 8156 8306
## goby_number2 7773 9380
## shrimp_speciesmannarensis 5688 7345
## shrimp_speciessciolii 6744 7971
## shrimp_speciesunknown 5985 6582
## stroke4:stagePost 8859 9867
## strokecontrol:stagePost 8952 9915
## stroke4:stagePre 8574 9542
## strokecontrol:stagePre 7872 9605
## phi_stroke4 5297 7681
## phi_strokecontrol 5756 7879
## phi_stagePost 6296 8223
## phi_stagePre 5761 7448
## phi_shrimp_number2 12351 9566
## phi_shrimp_number3 9521 9484
## phi_goby_number2 11371 9212
## phi_shrimp_speciesmannarensis 6299 8200
## phi_shrimp_speciessciolii 8625 9289
## phi_shrimp_speciesunknown 6123 7975
## phi_stroke4:stagePost 6525 8567
## phi_strokecontrol:stagePost 7024 9365
## phi_stroke4:stagePre 5675 8154
## phi_strokecontrol:stagePre 6562 8760
## zoi_stroke4 9704 8925
## zoi_strokecontrol 11802 9423
## zoi_stagePost 13321 10398
## zoi_stagePre 12967 9535
## zoi_shrimp_number2 8275 7292
## zoi_shrimp_number3 14852 10115
## zoi_goby_number2 11492 9895
## zoi_shrimp_speciesmannarensis 10434 8901
## zoi_shrimp_speciessciolii 11869 9910
## zoi_shrimp_speciesunknown 10623 9022
## zoi_stroke4:stagePost 12509 10233
## zoi_strokecontrol:stagePost 14058 9657
## zoi_stroke4:stagePre 12710 9427
## zoi_strokecontrol:stagePre 14262 10019
## coi_stroke4 19173 9416
## coi_strokecontrol 25479 9678
## coi_stagePost 21202 9999
## coi_stagePre 18921 9438
## coi_shrimp_number2 19814 9356
## coi_shrimp_number3 27342 8037
## coi_goby_number2 20652 9718
## coi_shrimp_speciesmannarensis 19191 9481
## coi_shrimp_speciessciolii 25735 8216
## coi_shrimp_speciesunknown 21686 8915
## coi_stroke4:stagePost 25821 9399
## coi_strokecontrol:stagePost 28409 8840
## coi_stroke4:stagePre 24599 9850
## coi_strokecontrol:stagePre 27646 7835
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
###2.4.1 Model fit and diagnostics ####2.4.1.1 Model convergence
extract_rhats(shrimp_out_model)
## Rhat
## Intercept 1.0009578
## phi_Intercept 1.0001122
## zoi_Intercept 0.9998407
## coi_Intercept 1.0003119
## stroke4 1.0005122
## strokecontrol 1.0005168
## stagePost 1.0010515
## stagePre 1.0002930
## shrimp_number2 0.9998931
## shrimp_number3 1.0001058
## goby_number2 1.0000986
## shrimp_speciesmannarensis 1.0008461
## shrimp_speciessciolii 1.0001949
## shrimp_speciesunknown 1.0006569
## stroke4:stagePost 1.0009191
## strokecontrol:stagePost 1.0007629
## stroke4:stagePre 1.0000436
## strokecontrol:stagePre 1.0003169
## phi_stroke4 1.0001378
## phi_strokecontrol 1.0003434
## phi_stagePost 1.0000806
## phi_stagePre 1.0001811
## phi_shrimp_number2 1.0000424
## phi_shrimp_number3 1.0004788
## phi_goby_number2 0.9999790
## phi_shrimp_speciesmannarensis 1.0008483
## phi_shrimp_speciessciolii 1.0000916
## phi_shrimp_speciesunknown 1.0001554
## phi_stroke4:stagePost 1.0000798
## phi_strokecontrol:stagePost 0.9999300
## phi_stroke4:stagePre 1.0000519
## phi_strokecontrol:stagePre 1.0001820
## zoi_stroke4 1.0000602
## zoi_strokecontrol 1.0000852
## zoi_stagePost 0.9999112
## zoi_stagePre 1.0004444
## zoi_shrimp_number2 0.9998895
## zoi_shrimp_number3 1.0002128
## zoi_goby_number2 0.9999529
## zoi_shrimp_speciesmannarensis 1.0000569
## zoi_shrimp_speciessciolii 1.0000006
## zoi_shrimp_speciesunknown 1.0000500
## zoi_stroke4:stagePost 1.0005978
## zoi_strokecontrol:stagePost 1.0003469
## zoi_stroke4:stagePre 1.0000262
## zoi_strokecontrol:stagePre 1.0000477
## coi_stroke4 1.0000408
## coi_strokecontrol 1.0003090
## coi_stagePost 1.0008316
## coi_stagePre 1.0007659
## coi_shrimp_number2 1.0003090
## coi_shrimp_number3 1.0009424
## coi_goby_number2 1.0005799
## coi_shrimp_speciesmannarensis 1.0006388
## coi_shrimp_speciessciolii 1.0003057
## coi_shrimp_speciesunknown 1.0000471
## coi_stroke4:stagePost 1.0001877
## coi_strokecontrol:stagePost 1.0010580
## coi_stroke4:stagePre 0.9999489
## coi_strokecontrol:stagePre 1.0011478
## site.sd(Intercept) 1.0005361
## site.sd(phi_Intercept) 1.0003269
## site.sd(zoi_Intercept) 1.0003958
## site.sd(coi_Intercept) 1.0002228
## site:burrow.sd(Intercept) 1.0007982
## site:burrow.sd(phi_Intercept) 1.0011935
## site:burrow.sd(zoi_Intercept) 1.0003181
## site:burrow.sd(coi_Intercept) 1.0001578
####2.4.1.2 Model fit
#posterior predictive checks
pp_check(shrimp_out_model)
pp_check(shrimp_out_model, type='stat', stat='mean')
#fitted vs predicted values
summary(fitted(shrimp_out_model))
## Estimate Est.Error Q2.5 Q97.5
## Min. :0.05903 Min. :0.04699 Min. :0.001233 Min. :0.2260
## 1st Qu.:0.19328 1st Qu.:0.08172 1st Qu.:0.070745 1st Qu.:0.3942
## Median :0.30089 Median :0.09558 Median :0.124544 Median :0.5096
## Mean :0.31422 Mean :0.09622 Mean :0.145423 Mean :0.5181
## 3rd Qu.:0.41508 3rd Qu.:0.10873 3rd Qu.:0.203817 3rd Qu.:0.6305
## Max. :0.70428 Max. :0.18746 Max. :0.454831 Max. :0.8784
summary(predict(shrimp_out_model))
## Estimate Est.Error Q2.5 Q97.5
## Min. :0.05957 Min. :0.1073 Min. :0.0000000 Min. :0.3900
## 1st Qu.:0.19325 1st Qu.:0.1952 1st Qu.:0.0000000 1st Qu.:0.7223
## Median :0.30309 Median :0.2252 Median :0.0000000 Median :0.8399
## Mean :0.31421 Mean :0.2262 Mean :0.0007245 Mean :0.8157
## 3rd Qu.:0.41371 3rd Qu.:0.2541 3rd Qu.:0.0000000 3rd Qu.:0.9288
## Max. :0.70185 Max. :0.4021 Max. :0.0307208 Max. :1.0000
###2.4.2 Posterior marginal effects ####2.4.2.1 Pairwise contrasts of stroke and stage
shrimp_out_contrasts <- shrimp_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
# Extract the relevant contrast information
shrimp_out_contrasts <- as.data.frame(shrimp_out_contrasts) |>
dplyr::filter(contrast %in% c("control Pre - 2 Pre", "control Pre - 4 Pre","4 Pre - 2 Pre", "control Pre - control During", "control Pre - control Post", "4 Pre - 4 During", "4 Pre - 4 Post", "2 Pre - 2 During", "2 Pre - 2 Post"))
shrimp_out_contrasts
## contrast estimate lower.HPD upper.HPD
## 2 Pre - 2 During 0.11124122 0.02991038 0.19526550
## 2 Pre - 2 Post 0.08161274 -0.00327231 0.16563823
## 4 Pre - 4 During 0.07146775 -0.03455554 0.18276195
## 4 Pre - 4 Post 0.00424298 -0.10889323 0.11883239
## 4 Pre - 2 Pre 0.00831290 -0.10911805 0.14158252
## control Pre - control During 0.02438722 -0.06822417 0.11723333
## control Pre - control Post -0.02058155 -0.11291257 0.06807063
## control Pre - 2 Pre -0.01089037 -0.13110502 0.11234895
## control Pre - 4 Pre -0.01896972 -0.15203909 0.11329580
##
## Results are averaged over the levels of: shrimp_number, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
shrimp_out_contrast_draws <- shrimp_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
re_formula = NA) |>
contrast(method = "revpairwise") |>
gather_emmeans_draws()
shrimp_out_temporal_contrasts <- shrimp_out_contrast_draws |>
dplyr::filter(contrast %in% c("control Pre - control During", "control Pre - control Post", "4 Pre - 4 During", "4 Pre - 4 Post", "2 Pre - 2 During", "2 Pre - 2 Post")) |>
mutate(stroke = ifelse(contrast %in% c("control Pre - control During", "control Pre - control Post"), "Control", ifelse(contrast %in% c("4 Pre - 4 During", "4 Pre - 4 Post"), "4", "2")))
shrimp_out_temporal_contrasts$stroke <- factor(shrimp_out_temporal_contrasts$stroke, levels = c("Control", "4", "2"))
shrimp_out_temporal_contrasts_fig <- ggplot(shrimp_out_temporal_contrasts, aes(x = .value, y = factor(contrast, levels = c("control Pre - control During", "control Pre - control Post", "4 Pre - 4 During", "4 Pre - 4 Post", "2 Pre - 2 During", "2 Pre - 2 Post")), color = stroke, fill = stroke)) +
geom_vline(xintercept = 0, color = "black", size = 0.5) +
stat_halfeye(.width = c(0.89, 0.95)) +
scale_color_manual(values = c(colour11, colour12, colour13)) +
scale_fill_manual(values = alpha(c(colour1, colour2, colour3), 0.8)) +
labs(x = "Difference in the proportion of time spent out of the burrow",
y = "Contrasts of the time periods relative to noise exposure",
fill = "Noise treatment", color = "Noise treatment",
fill_ramp = "Credible interval") +
coord_flip() +
theme(legend.position = "top",
axis.text.x = element_text(angle = 50, hjust =1)) +
scale_y_discrete(labels = c("Pre - During", "Pre - Post", "Pre - During", "Pre - Post", "Pre - During", "Pre - Post")) +
scale_x_continuous(breaks = seq(-0.2, 0.3, 0.1), labels = c("-0.2","-0.1" , "0", "0.1", "0.2", "0.3"), limits = c(-0.3, 0.4))
shrimp_out_temporal_contrasts_fig
###2.4.3 Extracting and plotting model estimates ####2.4.3.1 Extracting estimates
shrimp_out_model_estimates <- shrimp_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA)
shrimp_out_model_estimates <- as.data.frame(shrimp_out_model_estimates)
shrimp_out_model_estimates
## stroke stage emmean lower.HPD upper.HPD
## 2 During 0.2088068 0.1219075 0.3118424
## 4 During 0.2589222 0.1546599 0.3916332
## control During 0.2859898 0.1797559 0.4114104
## 2 Post 0.2386264 0.1385769 0.3589095
## 4 Post 0.3262621 0.1927346 0.4706111
## control Post 0.3334816 0.2124747 0.4633444
## 2 Pre 0.3217096 0.1933117 0.4538872
## 4 Pre 0.3306229 0.1925348 0.4863961
## control Pre 0.3110556 0.1939448 0.4440628
##
## Results are averaged over the levels of: shrimp_number, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
#Storing the model estimate means for each group to plot later on
shrimp_out_model_estimates_control <- shrimp_out_model_estimates |>
dplyr::filter(stroke == "control") |>
rename(.value = emmean)
shrimp_out_model_estimates_four <- shrimp_out_model_estimates |>
dplyr::filter(stroke == "4")|>
rename(.value = emmean)
shrimp_out_model_estimates_two <- shrimp_out_model_estimates |>
dplyr::filter(stroke == "2")|>
rename(.value = emmean)
####2.4.3.2 Plotting model estimates and raw data
shrimp_out_model_draws <- shrimp_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA) |>
gather_emmeans_draws()
shrimp_out_model_draws <- shrimp_out_model_draws |>
rename(emmean = .value)
shrimp_out_raw_data <- final_df |>
rename(emmean = shrimp_out)
shrimp_out_raw_data$stage <- factor(shrimp_out_raw_data$stage, levels = c("Pre", "During", "Post"))
shrimp_out_raw_data$stroke <- factor(shrimp_out_raw_data$stroke, levels = c("control", "4", "2"))
shrimp_out_model_estimates$stage <- factor(shrimp_out_model_estimates$stage, levels = c("Pre", "During", "Post"))
shrimp_out_model_estimates$stroke <- factor(shrimp_out_model_estimates$stroke, levels = c("control", "4", "2"))
shrimp_out_model_draws$stage <- factor(shrimp_out_model_draws$stage, levels = c("Pre", "During", "Post"))
shrimp_out_model_draws$stroke <- factor(shrimp_out_model_draws$stroke, levels = c("control", "4", "2"))
shrimp_out_means <- shrimp_out_raw_data |>
group_by(stage, stroke) |>
summarise(emmean = mean(emmean))
shrimp_out_fig_together <- ggplot(data = shrimp_out_raw_data, aes(x = stage, y = emmean, color = stroke, fill = stroke)) +
geom_point(data = shrimp_out_means, alpha = 1, position = position_dodgenudge(width = 0.8, x = -0.09), size = 9, shape = "_") +
geom_half_point(alpha = 0.3, position = position_dodgenudge(width = 0.8, x = -0.03), side = "l", range_scale = 0.5) +
geom_half_violin(data = shrimp_out_model_draws, aes(x = stage, y = emmean, color = stroke, fill = stroke), alpha = 0.6, position = position_dodge(width = 0.8), side = "r", size = 0) +
geom_point(data = shrimp_out_model_estimates, aes(x = stage, y = emmean, color = stroke), position = position_dodge(width = 0.8), size = 3) +
geom_errorbar(data = shrimp_out_model_estimates, aes(x = stage, y = emmean, ymin = lower.HPD, ymax = upper.HPD, color = stroke), position = position_dodge(width = 0.8), width = 0, size = 1.25) +
scale_fill_manual(values = c("control" = "#76A34A", "4" = "#F29E00", "2" = "#E25A00"), name = "Noise exposure", guide = guide_legend(override.aes = list(size = 3))) +
scale_color_manual(values = c("control" = "#2A330E", "4" = "#805213", "2" = "#542709"), name = "Noise exposure", guide = guide_legend(override.aes = list(size = 3))) +
scale_x_discrete(labels = c("Pre", "During", "Post")) +
labs(title = "The effects of boat noise on shrimp burrow use",
y = "Proportion of time in out of the burrow",
x = "Experimental phase relative to noise exposure") +
theme_classic() +
theme(
axis.text.x = element_text(size = 10, margin = margin(t = 5)),
axis.text.y = element_text(size = 10, margin = margin(r = 5)),
axis.title.x = element_text(size = 13, margin = margin(t = 10)),
axis.title.y = element_text(size = 13, margin = margin(r = 10)),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
legend.position = "right",
legend.justification = "top",
plot.title = element_text(hjust = 0.5, size = 16),
plot.margin = margin(10, 10, 10, 10)
)
shrimp_out_fig_together
####2.4.3.3 Plotting control estimates and raw data
shrimp_out_model_draws <- shrimp_out_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA) |>
gather_emmeans_draws()
shrimp_out_raw_data <- final_df |>
rename(.value = shrimp_out)
shrimp_out_model_draws$stage <- factor(shrimp_out_model_draws$stage, levels = c("Pre", "During", "Post"))
shrimp_out_model_draws$stroke <- factor(shrimp_out_model_draws$stroke, levels = c("control", "4", "2"))
shrimp_out_control_model_draws <- shrimp_out_model_draws |>
dplyr::filter(stroke == "control")
shrimp_control_raw_data <- final_df |>
dplyr::filter(stroke =="control") |>
rename(.value = shrimp_out)
shrimp_control_means <- shrimp_control_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
shrimp_control_raw_data$stage <- factor(shrimp_control_raw_data$stage, levels = c("Pre", "During", "Post"))
#assigning a value for the baseline mean
shrimp_out_model_estimates_control_intercept <- shrimp_out_model_estimates_control %>%
dplyr::filter(stage == "Pre") %>%
pull(.value)
shrimp_out_fig_control <- ggplot(shrimp_control_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = shrimp_out_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = shrimp_out_model_estimates_control, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = shrimp_out_model_estimates_control, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = shrimp_control_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
labs(x="Period relative to treatment",
y="Proportion of time out of burrow") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_text(size = 13),
axis.text.y = element_text(size = 13),
panel.grid = element_blank(),
axis.title.y = element_text(size = 16, face = "bold", margin = margin(t = 0, b = 0, l = 0, r = 10)),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#A3C088", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("Control")+
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#76A34A", "Post" = "#D3D9A7")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#2A330E", "Post" = "#5d634c"))
shrimp_out_fig_control
####2.4.3.4 Plotting 4-stroke estimates and raw data
shrimp_out_control_model_draws <- shrimp_out_model_draws |>
dplyr::filter(stroke == "4")
shrimp_four_raw_data <- final_df |>
dplyr::filter(stroke =="4") |>
rename(.value = shrimp_out)
shrimp_four_raw_data$stage <- factor(shrimp_four_raw_data$stage, levels = c("Pre", "During", "Post"))
shrimp_four_means <- shrimp_four_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
shrimp_out_fig_four <- ggplot(shrimp_four_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = shrimp_out_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = shrimp_out_model_estimates_four, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = shrimp_out_model_estimates_four, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = shrimp_four_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
xlab("Experimental phase relative to noise exposure")+
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.x = element_text(size = 13),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid = element_blank(),
axis.title.x = element_text(size = 16, face = "bold", margin = margin(t = 15, b = 0, l = 0, r = 0)),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#F6BC65", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("4-stroke")+
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#F29E00", "Post" = "#F2D5A0")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#805213", "Post" = "#70532b"))
shrimp_out_fig_four
####2.4.3.5 Plotting 2-stroke estimates and raw data
shrimp_out_control_model_draws <- shrimp_out_model_draws |>
dplyr::filter(stroke == "2")
shrimp_two_raw_data <- final_df |>
dplyr::filter(stroke =="2") |>
rename(.value = shrimp_out)
shrimp_two_raw_data$stage <- factor(shrimp_two_raw_data$stage, levels = c("Pre", "During", "Post"))
shrimp_two_means <- shrimp_two_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
shrimp_out_fig_two <- ggplot(shrimp_two_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = shrimp_out_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = shrimp_out_model_estimates_two, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = shrimp_out_model_estimates_two, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = shrimp_two_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_text(size = 13),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid = element_blank(),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#EB9063", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("2-stroke") +
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#E25A00", "Post" = "#F2B389")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#542709", "Post" = "#66422a"))
shrimp_out_fig_two
####2.4.3.6 Combining plots into one figure
shrimp_out_plot3 <- shrimp_out_fig_control + shrimp_out_fig_four + shrimp_out_fig_two +
plot_layout(ncol = 3)
shrimp_out_plot3
ggsave("./3-figs/shrimp_out_plot3.jpeg", plot = shrimp_out_plot3, units = "px", width = 2200, height = 1400, dpi = 300)
###2.4.4 Effects of covariates on the proportion of time shrimp spent outside of the burrow ####2.4.4.1 Goby number
shrimp_out_goby_number_contrasts <- shrimp_out_model |>
emmeans(~ goby_number,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
shrimp_out_goby_number_contrasts
## contrast estimate lower.HPD upper.HPD
## goby_number2 - goby_number1 0.0667 -0.019 0.16
##
## Results are averaged over the levels of: stroke, stage, shrimp_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
####2.4.4.2 Shrimp number
shrimp_out_shrimp_number_contrasts <- shrimp_out_model |>
emmeans(~ shrimp_number,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
shrimp_out_shrimp_number_contrasts <- as.data.frame(shrimp_out_shrimp_number_contrasts)
shrimp_out_shrimp_number_estimates <- shrimp_out_model |>
emmeans(~ shrimp_number,
epred = TRUE,
re_formla = NULL,
rg.limit =14000)
shrimp_out_shrimp_number_estimates
## shrimp_number emmean lower.HPD upper.HPD
## 1 0.223 0.129 0.329
## 2 0.339 0.262 0.433
## 3 0.308 0.124 0.534
##
## Results are averaged over the levels of: stroke, stage, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
shrimp_out_shrimp_number_contrast_draws <- shrimp_out_model |>
emmeans(~ shrimp_number,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise") |>
gather_emmeans_draws()
ggplot(data = shrimp_out_shrimp_number_contrast_draws, aes(x = contrast, y = .value )) +
stat_halfeye()
####2.4.4.3 Shrimp species
shrimp_out_shrimp_species_contrasts <- shrimp_out_model |>
emmeans(~ shrimp_species,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
shrimp_out_shrimp_species_contrasts
## contrast estimate lower.HPD upper.HPD
## mannarensis - bellulus 0.00168 -0.1752 0.1547
## sciolii - bellulus -0.09937 -0.3017 0.0899
## sciolii - mannarensis -0.10302 -0.2208 0.0390
## unknown - bellulus -0.01135 -0.2049 0.1599
## unknown - mannarensis -0.01330 -0.1033 0.0914
## unknown - sciolii 0.08881 -0.0588 0.2365
##
## Results are averaged over the levels of: stroke, stage, shrimp_number, goby_number
## Point estimate displayed: median
## HPD interval probability: 0.95
shrimp_out_shrimp_species_contrasts_draws <- shrimp_out_model |>
emmeans(~ shrimp_species,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise") |>
gather_emmeans_draws()
ggplot(data = shrimp_out_shrimp_species_contrasts_draws, aes(x = contrast, y = .value )) +
stat_halfeye()
##2.5 Modelling the proportion of time shrimp spent in contact with the goby ###2.5.1 Checking the response variable distribution
shrimp_contact_df <- final_df |>
mutate(is_zero = shrimp_contact == 0)
shrimp_contact_plot <- shrimp_contact_df |>
mutate(shrimp_contact = ifelse(is_zero, -0.1, shrimp_contact)) |>
ggplot(aes(x = shrimp_contact)) +
geom_histogram(aes(fill = is_zero), binwidth = 0.05,
boundary = 0, color = "white") +
geom_vline(xintercept = 0)
shrimp_contact_plot
# Distribution is a zero-one inflated beta distribution, because the value for shrimp_contact is bounded by 0-1 as it is a proportion of the total time and there is values of 0 and 1.
###2.5.2 Running the model
shrimp_contact_model <- brm(
bf(shrimp_contact ~ stroke*stage + shrimp_number + goby_number + shrimp_species + (1|site/burrow),
phi ~ stroke*stage + shrimp_number + goby_number + shrimp_species + (1|site/burrow),
zoi ~ stroke*stage + shrimp_number + goby_number + shrimp_species + (1|site/burrow),
coi ~ stroke*stage + shrimp_number + goby_number + shrimp_species + (1|site/burrow)
),
data = final_df,
prior = priors,
family = zero_one_inflated_beta(),
control = list(adapt_delta = 0.99, max_treedepth = 12),
chains = 4, iter = 6000, warmup = 1000, seed = 1, backend = "cmdstanr", cores = 4, threads = threading(2),
file = "shrimp_contact_model"
)
summary(shrimp_contact_model)
## Family: zero_one_inflated_beta
## Links: mu = logit; phi = log; zoi = logit; coi = logit
## Formula: shrimp_contact ~ stroke * stage + shrimp_number + goby_number + shrimp_species + (1 | site/burrow)
## phi ~ stroke * stage + shrimp_number + goby_number + shrimp_species + (1 | site/burrow)
## zoi ~ stroke * stage + shrimp_number + goby_number + shrimp_species + (1 | site/burrow)
## coi ~ stroke * stage + shrimp_number + goby_number + shrimp_species + (1 | site/burrow)
## Data: final_df (Number of observations: 369)
## Draws: 4 chains, each with iter = 6000; warmup = 1000; thin = 1;
## total post-warmup draws = 20000
##
## Multilevel Hyperparameters:
## ~site (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.31 0.28 0.01 1.01 1.00 4422 6853
## sd(phi_Intercept) 0.29 0.28 0.01 0.99 1.00 5122 6538
## sd(zoi_Intercept) 1.20 0.86 0.08 3.34 1.00 5704 6800
## sd(coi_Intercept) 2.44 1.97 0.09 7.26 1.00 7853 7421
##
## ~site:burrow (Number of levels: 123)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.81 0.10 0.63 1.01 1.00 4557 7393
## sd(phi_Intercept) 0.70 0.17 0.37 1.03 1.00 3649 4764
## sd(zoi_Intercept) 2.66 0.61 1.64 4.04 1.00 6766 11703
## sd(coi_Intercept) 7.19 3.93 2.17 16.86 1.00 8258 8479
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 0.79 0.55 -0.25 1.89 1.00
## phi_Intercept 1.12 0.67 -0.19 2.47 1.00
## zoi_Intercept -4.78 1.94 -8.72 -1.11 1.00
## coi_Intercept 1.38 3.94 -5.97 9.56 1.00
## stroke4 0.65 0.29 0.09 1.22 1.00
## strokecontrol 0.29 0.28 -0.26 0.83 1.00
## stagePost 0.26 0.18 -0.09 0.62 1.00
## stagePre 0.12 0.17 -0.21 0.46 1.00
## shrimp_number2 -0.29 0.26 -0.79 0.21 1.00
## shrimp_number3 -0.44 0.59 -1.61 0.73 1.00
## goby_number2 0.09 0.25 -0.40 0.57 1.00
## shrimp_speciesmannarensis 0.06 0.46 -0.86 0.97 1.00
## shrimp_speciessciolii -0.74 0.58 -1.89 0.37 1.00
## shrimp_speciesunknown -0.27 0.51 -1.28 0.73 1.00
## stroke4:stagePost 0.00 0.27 -0.53 0.54 1.00
## strokecontrol:stagePost -0.07 0.28 -0.62 0.46 1.00
## stroke4:stagePre -0.27 0.30 -0.86 0.32 1.00
## strokecontrol:stagePre 0.12 0.29 -0.45 0.68 1.00
## phi_stroke4 -0.63 0.44 -1.49 0.21 1.00
## phi_strokecontrol -0.56 0.44 -1.43 0.31 1.00
## phi_stagePost -0.46 0.39 -1.24 0.29 1.00
## phi_stagePre -0.16 0.41 -0.97 0.64 1.00
## phi_shrimp_number2 0.63 0.31 0.03 1.24 1.00
## phi_shrimp_number3 1.21 0.83 -0.48 2.78 1.00
## phi_goby_number2 -0.08 0.31 -0.69 0.52 1.00
## phi_shrimp_speciesmannarensis 0.51 0.53 -0.55 1.55 1.00
## phi_shrimp_speciessciolii 0.05 0.68 -1.27 1.38 1.00
## phi_shrimp_speciesunknown 0.56 0.59 -0.59 1.72 1.00
## phi_stroke4:stagePost 1.91 0.66 0.69 3.29 1.00
## phi_strokecontrol:stagePost 0.81 0.59 -0.32 1.96 1.00
## phi_stroke4:stagePre 0.45 0.59 -0.71 1.59 1.00
## phi_strokecontrol:stagePre 0.01 0.58 -1.14 1.13 1.00
## zoi_stroke4 0.68 1.08 -1.43 2.85 1.00
## zoi_strokecontrol 0.68 1.08 -1.45 2.83 1.00
## zoi_stagePost 0.33 0.77 -1.17 1.84 1.00
## zoi_stagePre -0.37 0.80 -1.97 1.17 1.00
## zoi_shrimp_number2 -1.12 0.86 -2.86 0.54 1.00
## zoi_shrimp_number3 -2.46 2.25 -7.16 1.75 1.00
## zoi_goby_number2 -1.17 1.00 -3.27 0.70 1.00
## zoi_shrimp_speciesmannarensis 1.58 1.59 -1.44 4.83 1.00
## zoi_shrimp_speciessciolii -0.21 1.98 -4.16 3.60 1.00
## zoi_shrimp_speciesunknown 1.58 1.66 -1.65 4.96 1.00
## zoi_stroke4:stagePost 0.05 1.08 -2.07 2.17 1.00
## zoi_strokecontrol:stagePost -0.31 1.12 -2.50 1.90 1.00
## zoi_stroke4:stagePre 0.72 1.11 -1.44 2.94 1.00
## zoi_strokecontrol:stagePre 0.76 1.13 -1.48 2.96 1.00
## coi_stroke4 0.83 2.44 -4.06 5.54 1.00
## coi_strokecontrol 1.57 2.48 -3.43 6.34 1.00
## coi_stagePost -0.24 1.76 -3.70 3.26 1.00
## coi_stagePre 2.13 2.02 -1.70 6.28 1.00
## coi_shrimp_number2 1.40 2.36 -3.32 5.98 1.00
## coi_shrimp_number3 -0.04 3.00 -5.98 5.83 1.00
## coi_goby_number2 -0.14 2.52 -5.10 4.81 1.00
## coi_shrimp_speciesmannarensis -0.13 2.57 -5.22 4.93 1.00
## coi_shrimp_speciessciolii 0.49 2.89 -5.19 6.19 1.00
## coi_shrimp_speciesunknown -0.27 2.59 -5.40 4.82 1.00
## coi_stroke4:stagePost 0.40 2.58 -4.65 5.54 1.00
## coi_strokecontrol:stagePost 1.88 2.36 -2.66 6.61 1.00
## coi_stroke4:stagePre -1.58 2.64 -6.70 3.77 1.00
## coi_strokecontrol:stagePre 1.38 2.68 -3.79 6.76 1.00
## Bulk_ESS Tail_ESS
## Intercept 6931 10771
## phi_Intercept 7413 10243
## zoi_Intercept 11633 13219
## coi_Intercept 14488 13507
## stroke4 7744 11935
## strokecontrol 7915 11314
## stagePost 10543 13593
## stagePre 11260 13581
## shrimp_number2 9739 12868
## shrimp_number3 10393 12786
## goby_number2 8642 12549
## shrimp_speciesmannarensis 6838 10512
## shrimp_speciessciolii 7494 11358
## shrimp_speciesunknown 6544 11073
## stroke4:stagePost 11045 13645
## strokecontrol:stagePost 10755 14071
## stroke4:stagePre 11336 14368
## strokecontrol:stagePre 11378 13322
## phi_stroke4 8909 12654
## phi_strokecontrol 8540 12610
## phi_stagePost 8783 12203
## phi_stagePre 7984 11110
## phi_shrimp_number2 13965 14698
## phi_shrimp_number3 13928 15287
## phi_goby_number2 13076 15075
## phi_shrimp_speciesmannarensis 8798 12296
## phi_shrimp_speciessciolii 10314 12900
## phi_shrimp_speciesunknown 8926 12135
## phi_stroke4:stagePost 7770 9350
## phi_strokecontrol:stagePost 9489 13376
## phi_stroke4:stagePre 9042 12575
## phi_strokecontrol:stagePre 8918 12709
## zoi_stroke4 11241 13071
## zoi_strokecontrol 11641 14480
## zoi_stagePost 14076 15104
## zoi_stagePre 14115 15298
## zoi_shrimp_number2 11327 12262
## zoi_shrimp_number3 22346 14674
## zoi_goby_number2 12953 13675
## zoi_shrimp_speciesmannarensis 13225 13825
## zoi_shrimp_speciessciolii 16521 15458
## zoi_shrimp_speciesunknown 13672 13882
## zoi_stroke4:stagePost 15603 15941
## zoi_strokecontrol:stagePost 15126 14384
## zoi_stroke4:stagePre 15084 16658
## zoi_strokecontrol:stagePre 15955 15657
## coi_stroke4 19934 16256
## coi_strokecontrol 19540 16407
## coi_stagePost 23374 15629
## coi_stagePre 20406 15606
## coi_shrimp_number2 18466 14794
## coi_shrimp_number3 36617 14014
## coi_goby_number2 24408 16540
## coi_shrimp_speciesmannarensis 24686 15542
## coi_shrimp_speciessciolii 34026 15571
## coi_shrimp_speciesunknown 24079 14807
## coi_stroke4:stagePost 28690 16315
## coi_strokecontrol:stagePost 26483 15080
## coi_stroke4:stagePre 23479 15804
## coi_strokecontrol:stagePre 25369 14685
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
###2.5.3 Model fit and diagnostics ####2.5.3.1 Model convergence
extract_rhats(shrimp_contact_model)
## Rhat
## Intercept 1.0005227
## phi_Intercept 1.0001435
## zoi_Intercept 0.9999331
## coi_Intercept 1.0000882
## stroke4 1.0004395
## strokecontrol 1.0003713
## stagePost 0.9999568
## stagePre 1.0002299
## shrimp_number2 1.0007670
## shrimp_number3 1.0004739
## goby_number2 1.0003800
## shrimp_speciesmannarensis 1.0003644
## shrimp_speciessciolii 1.0003243
## shrimp_speciesunknown 1.0003862
## stroke4:stagePost 0.9999707
## strokecontrol:stagePost 1.0000321
## stroke4:stagePre 1.0000330
## strokecontrol:stagePre 1.0000546
## phi_stroke4 1.0002337
## phi_strokecontrol 1.0001813
## phi_stagePost 1.0001259
## phi_stagePre 0.9998986
## phi_shrimp_number2 1.0002640
## phi_shrimp_number3 1.0003272
## phi_goby_number2 1.0001709
## phi_shrimp_speciesmannarensis 1.0003540
## phi_shrimp_speciessciolii 1.0001288
## phi_shrimp_speciesunknown 1.0001611
## phi_stroke4:stagePost 1.0001231
## phi_strokecontrol:stagePost 1.0001511
## phi_stroke4:stagePre 1.0000056
## phi_strokecontrol:stagePre 0.9999452
## zoi_stroke4 1.0003425
## zoi_strokecontrol 1.0000268
## zoi_stagePost 1.0002708
## zoi_stagePre 1.0003081
## zoi_shrimp_number2 1.0001342
## zoi_shrimp_number3 1.0002038
## zoi_goby_number2 1.0003639
## zoi_shrimp_speciesmannarensis 0.9999935
## zoi_shrimp_speciessciolii 1.0000000
## zoi_shrimp_speciesunknown 1.0000693
## zoi_stroke4:stagePost 1.0001350
## zoi_strokecontrol:stagePost 1.0003289
## zoi_stroke4:stagePre 1.0000980
## zoi_strokecontrol:stagePre 1.0001916
## coi_stroke4 1.0001010
## coi_strokecontrol 1.0001607
## coi_stagePost 1.0001902
## coi_stagePre 1.0001910
## coi_shrimp_number2 1.0001098
## coi_shrimp_number3 1.0005204
## coi_goby_number2 1.0002211
## coi_shrimp_speciesmannarensis 0.9999856
## coi_shrimp_speciessciolii 1.0002105
## coi_shrimp_speciesunknown 1.0001963
## coi_stroke4:stagePost 1.0004001
## coi_strokecontrol:stagePost 1.0001634
## coi_stroke4:stagePre 1.0001430
## coi_strokecontrol:stagePre 1.0000101
## site.sd(Intercept) 1.0004136
## site.sd(phi_Intercept) 1.0005560
## site.sd(zoi_Intercept) 1.0005406
## site.sd(coi_Intercept) 1.0001427
## site:burrow.sd(Intercept) 1.0003338
## site:burrow.sd(phi_Intercept) 1.0008470
## site:burrow.sd(zoi_Intercept) 1.0001670
## site:burrow.sd(coi_Intercept) 1.0002173
####2.5.3.2 Model fit
#posterior predictive checks
pp_check(shrimp_contact_model)
pp_check(shrimp_contact_model, type='stat', stat='mean')
#fitted vs predicted values
summary(fitted(shrimp_contact_model))
## Estimate Est.Error Q2.5 Q97.5
## Min. :0.2926 Min. :0.06588 Min. :0.03164 Min. :0.5580
## 1st Qu.:0.5898 1st Qu.:0.08469 1st Qu.:0.35528 1st Qu.:0.7970
## Median :0.7144 Median :0.09551 Median :0.48913 Median :0.8763
## Mean :0.6882 Mean :0.10260 Mean :0.46080 Mean :0.8580
## 3rd Qu.:0.8059 3rd Qu.:0.10986 3rd Qu.:0.58302 3rd Qu.:0.9377
## Max. :0.9398 Max. :0.22480 Max. :0.73064 Max. :0.9991
summary(predict(shrimp_contact_model))
## Estimate Est.Error Q2.5 Q97.5
## Min. :0.2921 Min. :0.1464 Min. :0.00000 Min. :0.8419
## 1st Qu.:0.5906 1st Qu.:0.1926 1st Qu.:0.01054 1st Qu.:1.0000
## Median :0.7147 Median :0.2166 Median :0.10769 Median :1.0000
## Mean :0.6881 Mean :0.2309 Mean :0.13260 Mean :0.9942
## 3rd Qu.:0.8069 3rd Qu.:0.2523 3rd Qu.:0.21747 3rd Qu.:1.0000
## Max. :0.9400 Max. :0.4517 Max. :0.53511 Max. :1.0000
###2.5.2 Posterior marginal effects ####2.5.2.1 Pairwise contrasts of stroke and stage
shrimp_contact_contrasts <- shrimp_contact_model |>
emmeans(~ stroke*stage,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
overall_mean <- shrimp_contact_model |>
emmeans(~ 1, # 1 means "no factors," i.e. the grand mean
epred = TRUE, # estimate on the response scale
re_formla = NULL,
rg.limit = 14000)
# View the summary with posterior mean and credible intervals
overall_mean
## 1 emmean lower.HPD upper.HPD
## overall 0.671 0.543 0.789
##
## Results are averaged over the levels of: stroke, stage, shrimp_number, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
# Extract the relevant contrast information
shrimp_contact_contrasts <- as.data.frame(shrimp_contact_contrasts) |>
dplyr::filter(contrast %in% c("control Pre - 2 Pre", "control Pre - 4 Pre","4 Pre - 2 Pre", "control Pre - control During", "control Pre - control Post", "4 Pre - 4 During", "4 Pre - 4 Post", "2 Pre - 2 During", "2 Pre - 2 Post"))
shrimp_contact_contrasts
## contrast estimate lower.HPD upper.HPD
## 2 Pre - 2 During 0.02860599 -0.04343278 0.10629264
## 2 Pre - 2 Post -0.02610253 -0.10486712 0.06071043
## 4 Pre - 4 During -0.02770567 -0.12792281 0.07167053
## 4 Pre - 4 Post -0.07128725 -0.16326137 0.01637126
## 4 Pre - 2 Pre 0.07507622 -0.04516392 0.19730855
## control Pre - control During 0.05412175 -0.03917928 0.15588585
## control Pre - control Post 0.01348306 -0.08026895 0.10607465
## control Pre - 2 Pre 0.09079625 -0.03019632 0.21535561
## control Pre - 4 Pre 0.01527249 -0.11322530 0.14688317
##
## Results are averaged over the levels of: shrimp_number, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
shrimp_contact_contrasts_draws <- shrimp_contact_model |>
emmeans(~ stroke*stage,
epred = TRUE,
re_formula = NA) |>
contrast(method = "revpairwise") |>
gather_emmeans_draws()
shrimp_contact_temporal_contrasts <- shrimp_contact_contrasts_draws |>
dplyr::filter(contrast %in% c("control Pre - control During", "control Pre - control Post","control Post - control During" , "4 Pre - 4 During", "4 Pre - 4 Post", "4 Post - 4 During", "2 Pre - 2 During", "2 Pre - 2 Post", "2 Post - 2 During")) |>
mutate(stroke = ifelse(contrast %in% c("control Pre - control During", "control Pre - control Post","control Post - control During"), "Control", ifelse(contrast %in% c("4 Pre - 4 During", "4 Pre - 4 Post", "4 Post - 4 During"), "4", "2")))
shrimp_contact_temporal_contrasts$stroke <- factor(shrimp_contact_temporal_contrasts$stroke, levels = c("Control", "4", "2"))
shrimp_contact_temporal_fig<- ggplot(shrimp_contact_temporal_contrasts, aes(x = .value, y = factor(contrast, levels = c("control Pre - control During", "control Pre - control Post","control Post - control During" , "4 Pre - 4 During", "4 Pre - 4 Post", "4 Post - 4 During", "2 Pre - 2 During", "2 Pre - 2 Post", "2 Post - 2 During")), color = stroke, fill = stroke)) +
geom_vline(xintercept = 0, color = "black", size = 0.5) +
stat_halfeye(.width = c(0.89, 0.95)) +
scale_color_manual(values = c(colour11, colour12, colour13)) +
scale_fill_manual(values = alpha(c(colour1, colour2, colour3), 0.8)) +
labs(x = "Difference in the proportion of time spent out of the burrow",
y = "Contrasts of the time periods relative to noise exposure",
fill = "Noise treatment", color = "Noise treatment",
fill_ramp = "Credible interval") +
coord_flip() +
theme(legend.position = "top",
axis.text.x = element_text(angle = 50, hjust =1)) +
scale_y_discrete(labels = c("Pre - During", "Pre - Post","Post - During" , "Pre - During", "Pre - Post", "Post - During", "Pre - During", "Pre - Post", "Post - During")) +
scale_x_continuous(breaks = seq(-0.2, 0.3, 0.1), labels = c("-0.2","-0.1" , "0", "0.1", "0.2", "0.3"), limits = c(-0.2, 0.3))
shrimp_contact_temporal_fig
###2.5.3 Extracting and plotting model estimates ####2.5.3.1 Extracting estimates
shrimp_contact_model_estimates <- shrimp_contact_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA)
shrimp_contact_model |>
emmeans(~ 1,
epred = TRUE,
rg.limit = 14000,
re_formula = NA)
## 1 emmean lower.HPD upper.HPD
## overall 0.671 0.543 0.789
##
## Results are averaged over the levels of: stroke, stage, shrimp_number, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
shrimp_contact_model_estimates <- as.data.frame(shrimp_contact_model_estimates)
shrimp_contact_model_estimates
## stroke stage emmean lower.HPD upper.HPD
## 2 During 0.5833632 0.4287101 0.7201552
## 4 During 0.7175983 0.5615712 0.8568383
## control During 0.6488437 0.4991426 0.7897400
## 2 Post 0.6397286 0.4872690 0.7778386
## 4 Post 0.7623786 0.6225847 0.8801669
## control Post 0.6905964 0.5467257 0.8199855
## 2 Pre 0.6127777 0.4593987 0.7542849
## 4 Pre 0.6893009 0.5230593 0.8280634
## control Pre 0.7050247 0.5548278 0.8339999
##
## Results are averaged over the levels of: shrimp_number, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
#Storing the model estimate means for each group to plot later on
shrimp_contact_model_estimates_control <- shrimp_contact_model_estimates |>
dplyr::filter(stroke == "control") |>
rename(.value = emmean)
shrimp_contact_model_estimates_four <- shrimp_contact_model_estimates |>
dplyr::filter(stroke == "4")|>
rename(.value = emmean)
shrimp_contact_model_estimates_two <- shrimp_contact_model_estimates |>
dplyr::filter(stroke == "2")|>
rename(.value = emmean)
####2.5.3.2 Plotting model estimates and raw data
shrimp_contact_model_draws <- shrimp_contact_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA) |>
gather_emmeans_draws()
shrimp_contact_model_draws <- shrimp_contact_model_draws |>
rename(emmean = .value)
shrimp_contact_raw_data <- final_df |>
rename(emmean = shrimp_contact)
shrimp_contact_raw_data$stage <- factor(shrimp_contact_raw_data$stage, levels = c("Pre", "During", "Post"))
shrimp_contact_raw_data$stroke <- factor(shrimp_contact_raw_data$stroke, levels = c("control", "4", "2"))
shrimp_contact_model_estimates$stage <- factor(shrimp_contact_model_estimates$stage, levels = c("Pre", "During", "Post"))
shrimp_contact_model_estimates$stroke <- factor(shrimp_contact_model_estimates$stroke, levels = c("control", "4", "2"))
shrimp_contact_model_draws$stage <- factor(shrimp_contact_model_draws$stage, levels = c("Pre", "During", "Post"))
shrimp_contact_model_draws$stroke <- factor(shrimp_contact_model_draws$stroke, levels = c("control", "4", "2"))
shrimp_contact_means <- shrimp_contact_raw_data |>
group_by(stage, stroke) |>
summarise(emmean = mean(emmean))
shrimp_contact_fig_together <- ggplot(data = shrimp_contact_raw_data, aes(x = stage, y = emmean, color = stroke, fill = stroke)) +
geom_point(data = shrimp_contact_means, alpha = 1, position = position_dodgenudge(width = 0.8, x = -0.09), size = 9, shape = "_") +
geom_half_point(alpha = 0.3, position = position_dodgenudge(width = 0.8, x = -0.03), side = "l", range_scale = 0.5) +
geom_half_violin(data = shrimp_contact_model_draws, aes(x = stage, y = emmean, color = stroke, fill = stroke), alpha = 0.6, position = position_dodge(width = 0.8), side = "r", size = 0) +
geom_point(data = shrimp_contact_model_estimates, aes(x = stage, y = emmean, color = stroke), position = position_dodge(width = 0.8), size = 3) +
geom_errorbar(data = shrimp_contact_model_estimates, aes(x = stage, y = emmean, ymin = lower.HPD, ymax = upper.HPD, color = stroke), position = position_dodge(width = 0.8), width = 0, size = 1.25) +
scale_fill_manual(values = c("control" = "#76A34A", "4" = "#F29E00", "2" = "#E25A00"), name = "Noise exposure", guide = guide_legend(override.aes = list(size = 3))) +
scale_color_manual(values = c("control" = "#2A330E", "4" = "#805213", "2" = "#542709"), name = "Noise exposure", guide = guide_legend(override.aes = list(size = 3))) +
scale_x_discrete(labels = c("Pre", "During", "Post")) +
labs(title = "The effects of boat noise on goby-shrimp contact",
y = "Proportion of time the in contact",
x = "Experimental phase relative to noise exposure") +
theme_classic() +
theme(
axis.text.x = element_text(size = 10, margin = margin(t = 5)),
axis.text.y = element_text(size = 10, margin = margin(r = 5)),
axis.title.x = element_text(size = 13, margin = margin(t = 10)),
axis.title.y = element_text(size = 13, margin = margin(r = 10)),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
legend.position = "right",
legend.justification = "top",
plot.title = element_text(hjust = 0.5, size = 16),
plot.margin = margin(10, 10, 10, 10)
)
shrimp_contact_fig_together
ggsave("./3-figs/shrimp_contact_fig_together.png", plot = shrimp_contact_fig_together, units = "px", width = 2200, height = 1400, dpi = 300)
####2.5.3.3 Plotting control estimates and raw data
shrimp_contact_model_draws <- shrimp_contact_model |>
emmeans(~ stroke*stage,
epred = TRUE,
rg.limit = 14000,
re_formula = NA) |>
gather_emmeans_draws()
shrimp_contact_raw_data <- final_df |>
rename(.value = shrimp_contact)
shrimp_contact_model_draws$stage <- factor(shrimp_contact_model_draws$stage, levels = c("Pre", "During", "Post"))
shrimp_contact_model_draws$stroke <- factor(shrimp_contact_model_draws$stroke, levels = c("control", "4", "2"))
shrimp_contact_control_model_draws <- shrimp_contact_model_draws |>
dplyr::filter(stroke == "control")
goby_control_raw_data <- final_df |>
dplyr::filter(stroke =="control") |>
rename(.value = shrimp_contact)
goby_control_means <- goby_control_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
goby_control_raw_data$stage <- factor(goby_control_raw_data$stage, levels = c("Pre", "During", "Post"))
#assigning a value for the baseline mean
shrimp_contact_model_estimates_control_intercept <- shrimp_contact_model_estimates_control %>%
dplyr::filter(stage == "Pre") %>%
pull(.value)
shrimp_contact_fig_control <- ggplot(goby_control_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = shrimp_contact_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = shrimp_contact_model_estimates_control, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = shrimp_contact_model_estimates_control, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = goby_control_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
labs(x="Period relative to treatment",
y="Proportion of time in contact") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_text(size = 13),
axis.text.y = element_text(size = 13),
axis.title.y = element_text(size = 16, face = "bold", margin = margin(t = 0, b = 0, l = 0, r = 10)),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#A3C088", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("Control")+
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#76A34A", "Post" = "#D3D9A7")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#2A330E", "Post" = "#5d634c"))
shrimp_contact_fig_control
####2.5.3.4 Plotting 4-stroke estimates and raw data
shrimp_contact_control_model_draws <- shrimp_contact_model_draws |>
dplyr::filter(stroke == "4")
goby_four_raw_data <- final_df |>
dplyr::filter(stroke =="4") |>
rename(.value = shrimp_contact)
goby_four_raw_data$stage <- factor(goby_four_raw_data$stage, levels = c("Pre", "During", "Post"))
goby_four_means <- goby_four_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
shrimp_contact_fig_four <- ggplot(goby_four_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = shrimp_contact_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = shrimp_contact_model_estimates_four, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = shrimp_contact_model_estimates_four, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = goby_four_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
xlab("Experimental phase relative to noise exposure")+
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.x = element_text(size = 13),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.x = element_text(size = 16, face = "bold", margin = margin(t = 15, b = 0, l = 0, r = 0)),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#F6BC65", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("4-stroke")+
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#F29E00", "Post" = "#F2D5A0")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#805213", "Post" = "#70532b"))
shrimp_contact_fig_four
####2.5.3.5 Plotting 2-stroke estimates and raw data
shrimp_contact_control_model_draws <- shrimp_contact_model_draws |>
dplyr::filter(stroke == "2")
goby_two_raw_data <- final_df |>
dplyr::filter(stroke =="2") |>
rename(.value = shrimp_contact)
goby_two_raw_data$stage <- factor(goby_two_raw_data$stage, levels = c("Pre", "During", "Post"))
goby_two_means <- goby_two_raw_data |>
group_by(stage) |>
summarise(.value = mean(.value))
shrimp_contact_fig_two <- ggplot(goby_two_raw_data, aes(y = .value, x = stage, fill = stage, color = stage)) +
geom_half_violin(data = shrimp_contact_control_model_draws, alpha = 0.7, size = 0, side = "r", position = position_nudge(x = 0.2)) +
geom_jitter(position = position_jitternudge(x = -0.1, width = 0.1, nudge.from = "jittered"), alpha = 0.4, size = 2) +
geom_crossbar(alpha = 0.6, size = 0.8, data = shrimp_contact_model_estimates_two, aes(ymin = lower.HPD, ymax = upper.HPD), width = 0, position = position_nudge(x = 0.2))+
geom_point(pch = 23, stroke = 1.2, data = shrimp_contact_model_estimates_two, size = 2, position = position_nudge(x = 0.2)) +
geom_point(data = goby_two_means, position = position_dodgenudge(width = 0.8, x = -0.1), size = 9, shape = "_") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_text(size = 13),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.title = element_textbox_simple(face = "bold", size = 16, halign = 0.5, linewidth = 0.3, linetype = 1, color = "black", fill = "#EB9063", padding = margin(5, 5, 5, 5)),
panel.background = element_rect(fill = "white", colour = "black"),
plot.background = element_rect(fill = NA, colour = NA))+
ggtitle("2-stroke") +
scale_fill_manual(values = c("Pre" = "darkgrey", "During" = "#E25A00", "Post" = "#F2B389")) +
scale_color_manual(values = c("Pre" = "#242424", "During" = "#542709", "Post" = "#66422a"))
shrimp_contact_fig_two
####2.5.3.6 Combining plots into one figure
shrimp_contact_plot3 <- shrimp_contact_fig_control + shrimp_contact_fig_four + shrimp_contact_fig_two +
plot_layout(ncol = 3)
shrimp_contact_plot3
ggsave("./3-figs/shrimp_contact_plot3.jpeg", plot = shrimp_contact_plot3, units = "px", width = 2200, height = 1400, dpi = 300)
###2.5.4 Modelling the effect of covariates on the proportion of time spent in contact ####2.5.4.1 Goby number
shrimp_contact_goby_number_contrasts <- shrimp_contact_model |>
emmeans(~ goby_number,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
shrimp_contact_goby_number_contrasts
## contrast estimate lower.HPD upper.HPD
## goby_number2 - goby_number1 0.0147 -0.0847 0.109
##
## Results are averaged over the levels of: stroke, stage, shrimp_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
####2.5.4.2 Shrimp number
shrimp_contact_shrimp_number_contrasts <- shrimp_contact_model |>
emmeans(~ shrimp_number,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
shrimp_contact_shrimp_number_contrasts
## contrast estimate lower.HPD upper.HPD
## shrimp_number2 - shrimp_number1 -0.0534 -0.152 0.0425
## shrimp_number3 - shrimp_number1 -0.0867 -0.336 0.1326
## shrimp_number3 - shrimp_number2 -0.0343 -0.270 0.1739
##
## Results are averaged over the levels of: stroke, stage, goby_number, shrimp_species
## Point estimate displayed: median
## HPD interval probability: 0.95
####2.5.4.2 Shrimp species
shrimp_contact_shrimp_species_contrasts <- shrimp_contact_model |>
emmeans(~ shrimp_species,
epred = TRUE,
re_formla = NULL,
rg.limit =14000) |>
contrast(method = "revpairwise")
shrimp_contact_shrimp_species_contrasts
## contrast estimate lower.HPD upper.HPD
## mannarensis - bellulus 0.0115 -0.147 0.1861
## sciolii - bellulus -0.1492 -0.377 0.0839
## sciolii - mannarensis -0.1635 -0.343 0.0107
## unknown - bellulus -0.0504 -0.235 0.1429
## unknown - mannarensis -0.0625 -0.170 0.0338
## unknown - sciolii 0.1002 -0.080 0.2967
##
## Results are averaged over the levels of: stroke, stage, shrimp_number, goby_number
## Point estimate displayed: median
## HPD interval probability: 0.95
#3 Manuscript output files ##3.1 Supp tables ###3.1.1 Model estimates
goby_out_model_estimates <- goby_out_model_estimates |>
mutate(`Behavioural measure` = "Proportion of time gobies spent outside of the burrow")
shrimp_out_model_estimates <- shrimp_out_model_estimates |>
mutate(`Behavioural measure` = "Proportion of time shrimp spent outside of the burrow")
shrimp_contact_model_estimates <- shrimp_contact_model_estimates |>
mutate(`Behavioural measure` = "Proportion of time shrimp spent in contact with a goby")
model_estimates <- bind_rows(goby_out_model_estimates,
shrimp_out_model_estimates,
shrimp_contact_model_estimates)
model_estimates <- model_estimates |>
mutate(stroke = ifelse(stroke == "4", "4-stroke", ifelse(stroke == "2", "2-stroke", "control")),
stroke = factor(stroke, levels = c("control", "4-stroke", "2-stroke")),
stage = factor(stage, levels = c("Pre", "During", "Post")),
`Behavioural measure` = factor(`Behavioural measure`, levels = c("Proportion of time gobies spent outside of the burrow", "Proportion of time shrimp spent outside of the burrow", "Proportion of time shrimp spent in contact with a goby"))) |>
arrange(`Behavioural measure`, stage, stroke) |>
rename(`Estimated marginal mean` = emmean,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Phase = stage,
`Noise exposure` = stroke) |>
select(c(`Behavioural measure`, Phase, `Noise exposure`, `Estimated marginal mean`, `Lower 95% CrI`, `Upper 95% CrI`))
write_csv(model_estimates, "./2-data/model_estimates.csv")
###3.1.2 Pairwise contrasts
#goby out
goby_out_contrasts <- goby_out_results.table |>
mutate(`Predictor variable` = "Treatment and phase interaction",
contrast = str_replace_all(contrast, "\\b2\\b", "2-Stroke"),
contrast = str_replace_all(contrast, "\\b4\\b", "4-Stroke"),
contrast = str_replace_all(contrast, "\\bcontrol\\b", "Control")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
goby_out_goby_number1 <- goby_out_goby_number |>
as.data.frame() |>
mutate(`Predictor variable` = "Number of gobies",
contrast = str_replace_all(contrast, "\\bgoby_number2\\b", "2 gobies"),
contrast = str_replace_all(contrast, "\\bgoby_number1\\b", "1 goby")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
goby_out_shrimp_number1 <- goby_out_shrimp_number |>
as.data.frame() |>
mutate(`Predictor variable` = "Number of shrimp",
contrast = str_replace_all(contrast, "\\bshrimp_number2\\b", "2 shrimp"),
contrast = str_replace_all(contrast, "\\bshrimp_number1\\b", "1 shrimp"),
contrast = str_replace_all(contrast, "\\bshrimp_number3\\b", "3 shrimp")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
goby_out_shrimp_species1 <- goby_out_shrimp_species |>
as.data.frame() |>
mutate(`Predictor variable` = "Shrimp species"
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
goby_out_contrasts_table <- rbind(goby_out_contrasts, goby_out_goby_number1, goby_out_shrimp_number1, goby_out_shrimp_species1) |>
mutate(
`Estimated marginal difference` = round(`Estimated marginal difference`, 3),
`Lower 95% CrI` = round(`Lower 95% CrI`, 3),
`Upper 95% CrI` = round(`Upper 95% CrI`, 3)
)
write_csv(goby_out_contrasts_table, "./2-data/goby_out_contrasts_table.csv")
#shrimp out
shrimp_out_contrasts1 <- shrimp_out_contrasts |>
mutate(`Predictor variable` = "Treatment and phase interaction",
contrast = str_replace_all(contrast, "\\b2\\b", "2-Stroke"),
contrast = str_replace_all(contrast, "\\b4\\b", "4-Stroke"),
contrast = str_replace_all(contrast, "\\bcontrol\\b", "Control")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
shrimp_out_goby_number_contrasts1 <- shrimp_out_goby_number_contrasts |>
as.data.frame() |>
mutate(`Predictor variable` = "Number of gobies",
contrast = str_replace_all(contrast, "\\bgoby_number2\\b", "2 gobies"),
contrast = str_replace_all(contrast, "\\bgoby_number1\\b", "1 goby")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
shrimp_out_shrimp_number_contrasts1 <- shrimp_out_shrimp_number_contrasts |>
as.data.frame() |>
mutate(`Predictor variable` = "Number of shrimp",
contrast = str_replace_all(contrast, "\\bshrimp_number2\\b", "2 shrimp"),
contrast = str_replace_all(contrast, "\\bshrimp_number1\\b", "1 shrimp"),
contrast = str_replace_all(contrast, "\\bshrimp_number3\\b", "3 shrimp")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
shrimp_out_shrimp_species_contrasts1 <- shrimp_out_shrimp_species_contrasts |>
as.data.frame() |>
mutate(`Predictor variable` = "Shrimp species"
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
shrimp_out_contrasts_table <- rbind(shrimp_out_contrasts1, shrimp_out_goby_number_contrasts1, shrimp_out_shrimp_number_contrasts1, shrimp_out_shrimp_species_contrasts1) |>
mutate(
`Estimated marginal difference` = round(`Estimated marginal difference`, 3),
`Lower 95% CrI` = round(`Lower 95% CrI`, 3),
`Upper 95% CrI` = round(`Upper 95% CrI`, 3)
)
write_csv(shrimp_out_contrasts_table, "./2-data/shrimp_out_contrasts_table.csv")
#shrimp contact
shrimp_contact_contrasts1 <- shrimp_contact_contrasts |>
mutate(`Predictor variable` = "Treatment and phase interaction",
contrast = str_replace_all(contrast, "\\b2\\b", "2-Stroke"),
contrast = str_replace_all(contrast, "\\b4\\b", "4-Stroke"),
contrast = str_replace_all(contrast, "\\bcontrol\\b", "Control")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
shrimp_contact_goby_number_contrasts1 <- shrimp_contact_goby_number_contrasts |>
as.data.frame() |>
mutate(`Predictor variable` = "Number of gobies",
contrast = str_replace_all(contrast, "\\bgoby_number2\\b", "2 gobies"),
contrast = str_replace_all(contrast, "\\bgoby_number1\\b", "1 goby")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
shrimp_contact_shrimp_number_contrasts1 <- shrimp_contact_shrimp_number_contrasts |>
as.data.frame() |>
mutate(`Predictor variable` = "Number of shrimp",
contrast = str_replace_all(contrast, "\\bshrimp_number2\\b", "2 shrimp"),
contrast = str_replace_all(contrast, "\\bshrimp_number1\\b", "1 shrimp"),
contrast = str_replace_all(contrast, "\\bshrimp_number3\\b", "3 shrimp")
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
shrimp_contact_shrimp_species_contrasts1 <- shrimp_contact_shrimp_species_contrasts |>
as.data.frame() |>
mutate(`Predictor variable` = "Shrimp species"
) |>
rename(`Estimated marginal difference` = estimate,
`Lower 95% CrI` = lower.HPD,
`Upper 95% CrI` = upper.HPD,
Contrast = contrast) |>
select(c(`Predictor variable`, Contrast, `Estimated marginal difference`, `Lower 95% CrI`, `Upper 95% CrI`))
shrimp_contact_contrasts_table <- rbind(shrimp_contact_contrasts1, shrimp_contact_goby_number_contrasts1, shrimp_contact_shrimp_number_contrasts1, shrimp_contact_shrimp_species_contrasts1) |>
mutate(
`Estimated marginal difference` = round(`Estimated marginal difference`, 3),
`Lower 95% CrI` = round(`Lower 95% CrI`, 3),
`Upper 95% CrI` = round(`Upper 95% CrI`, 3)
)
write_csv(shrimp_contact_contrasts_table, "./2-data/shrimp_contact_contrasts_table.csv")
table_S1 <- final_df %>%
select(c("burrow", "stroke", "site")) |>
unique() |>
group_by(site, stroke) %>%
summarise(n = n(), .groups = "drop") %>%
pivot_wider(names_from = stroke, values_from = n, values_fill = 0) %>%
mutate(Total = `control` + `4` + `2`) %>%
arrange(site) %>%
rename(`4-stroke` = `4`,
`2-stroke` = `2`,
Control = control)